#include <define.h>

MODULE MOD_SnowSnicar

!-----------------------------------------------------------------------
! !DESCRIPTION:
!  Calculate albedo of snow containing impurities and the evolution of
!  snow effective radius
!
!  Original:
!  1) The Community Land Model version 5.0 (CLM5.0)
!  2) Energy Exascale Earth System Model version 2.0 (E3SM v2.0) Land
!     Model (ELM v2.0)
!
! !REFERENCES:
!  1) Flanner et al, 2021, SNICAR-ADv3: a community tool for modeling
!  spectral snow albedo.  Geosci. Model Dev., 14, 7673-7704,
!  https://doi.org/10.5194/gmd-14-7673-2021
!  2) Hao et al., 2023, Improving snow albedo modeling in the E3SM land
!  model (version 2.0) and assessing its impacts on snow and surface
!  fluxes over the Tibetan Plateau.  Geosci. Model Dev., 16, 75-94,
!  https://doi.org/10.5194/gmd-16-75-2023
!
! !REVISIONS:
!  Yongjiu Dai, and Hua Yuan, December, 2022 : ASSEMBLING and FITTING
!
!-----------------------------------------------------------------------
! !USES:
   USE MOD_Precision
   USE MOD_Vars_Global, only: maxsnl
   USE MOD_SPMD_Task

   IMPLICIT NONE
! SAVE
   real(R8),parameter :: SHR_CONST_PI     = 3.14159265358979323846_R8
   real(R8),parameter :: SHR_CONST_RHOICE = 0.917e3_R8           !  density of ice (kg/m^3)

   integer, parameter :: iulog = 6                               ! "stdout" log file unit number, default is 6
   integer, parameter :: numrad  = 2                             !  number of solar radiation bands: vis, nir

!--------------------------------------------------------------------
! DAI, Dec. 29, 2022
! Temporary setting

   logical, parameter :: use_extrasnowlayers = .false.
   character(len=256), parameter :: snow_shape = 'sphere'        ! (=1), 'spheroid'(=2), 'hexagonal_plate'(=3), 'koch_snowflake'(=4)
   logical, parameter :: use_dust_snow_internal_mixing = .false.
   character(len=256), parameter :: snicar_atm_type = 'default'  ! Atmospheric profile used to obtain surface-incident spectral flux distribution
                                                                 ! and subsequent broadband albedo
                                                                 ! = 'mid-latitude_winter'  ! => 1
                                                                 ! = 'mid-latitude_summer'  ! => 2
                                                                 ! = 'sub-Arctic_winter'    ! => 3
                                                                 ! = 'sub-Arctic_summer'    ! => 4
                                                                 ! = 'summit_Greenland'     ! => 5 (sub-Arctic summer, surface pressure of 796hPa)
                                                                 ! = 'high_mountain'        ! => 6 (summer, surface pressure of 556 hPa)
!DAI, Dec. 29, 2022
!-----------------------------------------------------------------------

   ! !PUBLIC MEMBER FUNCTIONS:
   PUBLIC :: SNICAR_RT        ! Snow albedo and vertically-resolved solar absorption
   PUBLIC :: SNICAR_AD_RT     ! Snow albedo and vertically-resolved solar absorption by adding-doubling solution
                              ! To USE this subroutine, set use_snicar_ad = true
   PUBLIC :: SnowAge_grain    ! Snow effective grain size evolution
   PUBLIC :: SnowAge_init     ! Initial read in of snow-aging file
   PUBLIC :: SnowOptics_init  ! Initial read in of snow-optics file
   !
   ! !PUBLIC DATA MEMBERS:
   integer,  PUBLIC, parameter :: sno_nbr_aer = 8         ! number of aerosol species in snowpack
                                                          ! (indices described above) [nbr]
   logical,  PUBLIC, parameter :: DO_SNO_OC   = .false.   ! parameter to include organic carbon (OC)
                                                          ! in snowpack radiative calculations
   logical,  PUBLIC, parameter :: DO_SNO_AER  = .true.    ! parameter to include aerosols in snowpack radiative calculations
   ! !PRIVATE DATA MEMBERS:
   integer,  parameter :: numrad_snw     = 5              ! number of spectral bands used in snow model [nbr]
   integer,  parameter :: nir_bnd_bgn    = 2              ! first band index in near-IR spectrum [idx]
   integer,  parameter :: nir_bnd_end    = 5              ! ending near-IR band index [idx]
   integer,  parameter :: idx_Mie_snw_mx = 1471           ! number of effective radius indices used in Mie lookup table [idx]
   integer,  parameter :: idx_T_max      = 11             ! maximum temperature index used in aging lookup table [idx]
   integer,  parameter :: idx_T_min      = 1              ! minimum temperature index used in aging lookup table [idx]
   integer,  parameter :: idx_Tgrd_max   = 31             ! maximum temperature gradient index used in aging lookup table [idx]
   integer,  parameter :: idx_Tgrd_min   = 1              ! minimum temperature gradient index used in aging lookup table [idx]
   integer,  parameter :: idx_rhos_max   = 8              ! maximum snow density index used in aging lookup table [idx]
   integer,  parameter :: idx_rhos_min   = 1              ! minimum snow density index used in aging lookup table [idx]

#ifdef MODAL_AER
   ! NOTE: right now the macro 'MODAL_AER' is not defined anywhere, i.e.,
   ! the below (modal aerosol scheme) is not available and can not be
   ! active either. It depends on the specific input aerosol deposition
   ! data which is suitable for modal scheme. [06/15/2023, Hua Yuan]
   !mgf++
   integer,  parameter :: idx_bc_nclrds_min     = 1       ! minimum index for BC particle size in optics lookup table
   integer,  parameter :: idx_bc_nclrds_max     = 10      ! maximum index for BC particle size in optics lookup table
   integer,  parameter :: idx_bcint_icerds_min  = 1       ! minimum index for snow grain size in optics lookup table for within-ice BC
   integer,  parameter :: idx_bcint_icerds_max  = 8       ! maximum index for snow grain size in optics lookup table for within-ice BC
   !mgf--
#endif

   integer,  parameter :: snw_rds_max_tbl = 1500          ! maximum effective radius defined in Mie lookup table [microns]
   integer,  parameter :: snw_rds_min_tbl = 30            ! minimum effective radius defined in Mie lookup table [microns]
   real(r8), parameter :: snw_rds_max     = 1500._r8      ! maximum allowed snow effective radius [microns]
   real(r8), parameter :: snw_rds_min     = 54.526_r8     ! minimum allowed snow effective radius (also "fresh snow" value) [microns
   real(r8), parameter :: snw_rds_refrz   = 1000._r8      ! effective radius of re-frozen snow [microns]
   real(r8), parameter :: min_snw = 1.0E-30_r8            ! minimum snow mass required for SNICAR RT calculation [kg m-2]
   !real(r8), parameter :: C1_liq_Brun89  = 1.28E-17_r8   ! constant for liquid water grain growth [m3 s-1],
                                                          ! from Brun89
   real(r8), parameter :: C1_liq_Brun89   = 0._r8         ! constant for liquid water grain growth [m3 s-1],
                                                          ! from Brun89: zeroed to accomodate dry snow aging
   real(r8), parameter :: C2_liq_Brun89   = 4.22E-13_r8   ! constant for liquid water grain growth [m3 s-1],
                                                          ! from Brun89: corrected for LWC in units of percent

   real(r8), parameter :: tim_cns_bc_rmv  = 2.2E-8_r8     ! time constant for removal of BC in snow on sea-ice
                                                          ! [s-1] (50% mass removal/year)
   real(r8), parameter :: tim_cns_oc_rmv  = 2.2E-8_r8     ! time constant for removal of OC in snow on sea-ice
                                                          ! [s-1] (50% mass removal/year)
   real(r8), parameter :: tim_cns_dst_rmv = 2.2E-8_r8     ! time constant for removal of dust in snow on sea-ice
                                                          ! [s-1] (50% mass removal/year)
   !$acc declare copyin(C1_liq_Brun89, C2_liq_Brun89, &
   !$acc tim_cns_bc_rmv, tim_cns_oc_rmv, tim_cns_dst_rmv)

   ! scaling of the snow aging rate (tuning option):
   logical :: flg_snoage_scl    = .false.                 ! flag for scaling the snow aging rate by some arbitrary factor
   real(r8), parameter :: xdrdt = 1.0_r8                  ! arbitrary factor applied to snow aging rate
   ! snow and aerosol Mie parameters:
   ! (arrays declared here, but are set in iniTimeConst)
   ! (idx_Mie_snw_mx is number of snow radii with defined parameters (i.e. from 30um to 1500um))

   ! direct-beam weighted ice optical properties
   real(r8), allocatable :: ss_alb_snw_drc     (:,:) ! (idx_Mie_snw_mx,numrad_snw);
   real(r8), allocatable :: asm_prm_snw_drc    (:,:) ! (idx_Mie_snw_mx,numrad_snw);
   real(r8), allocatable :: ext_cff_mss_snw_drc(:,:) ! (idx_Mie_snw_mx,numrad_snw);

   ! diffuse radiation weighted ice optical properties
   real(r8), allocatable :: ss_alb_snw_dfs     (:,:) ! (idx_Mie_snw_mx,numrad_snw);
   real(r8), allocatable :: asm_prm_snw_dfs    (:,:) ! (idx_Mie_snw_mx,numrad_snw);
   real(r8), allocatable :: ext_cff_mss_snw_dfs(:,:) ! (idx_Mie_snw_mx,numrad_snw);

   ! direct & diffuse flux
   real(r8), allocatable :: flx_wgt_dir (:,:,:) ! (6, 90, numrad_snw) ! direct flux, six atmospheric types, 0-89 SZA
   real(r8), allocatable :: flx_wgt_dif (:,:)   ! (6, numrad_snw)     ! diffuse flux, six atmospheric types

   ! snow grain shape
   integer, parameter :: snow_shape_sphere            = 1
   integer, parameter :: snow_shape_spheroid          = 2
   integer, parameter :: snow_shape_hexagonal_plate   = 3
   integer, parameter :: snow_shape_koch_snowflake    = 4

   ! atmospheric condition for SNICAR-AD
   integer, parameter :: atm_type_default             = 0
   integer, parameter :: atm_type_mid_latitude_winter = 1
   integer, parameter :: atm_type_mid_latitude_summer = 2
   integer, parameter :: atm_type_sub_Arctic_winter   = 3
   integer, parameter :: atm_type_sub_Arctic_summer   = 4
   integer, parameter :: atm_type_summit_Greenland    = 5
   integer, parameter :: atm_type_high_mountain       = 6

#ifdef MODAL_AER
   !mgf++
   ! Size-dependent BC optical properties. Currently a fixed BC size is
   ! assumed, but this framework enables optical properties to be
   ! assigned based on the BC effective radius, should this be
   ! implemented in the future.
   !
   ! within-ice BC (i.e., BC that was deposited within hydrometeors)
   real(r8), allocatable :: ss_alb_bc1     (:,:)  ! (numrad_snw,idx_bc_nclrds_max);
   real(r8), allocatable :: asm_prm_bc1    (:,:)  ! (numrad_snw,idx_bc_nclrds_max);
   real(r8), allocatable :: ext_cff_mss_bc1(:,:)  ! (numrad_snw,idx_bc_nclrds_max);

   ! external BC
   real(r8), allocatable :: ss_alb_bc2     (:,:)  ! (numrad_snw,idx_bc_nclrds_max);
   real(r8), allocatable :: asm_prm_bc2    (:,:)  ! (numrad_snw,idx_bc_nclrds_max);
   real(r8), allocatable :: ext_cff_mss_bc2(:,:)  ! (numrad_snw,idx_bc_nclrds_max);
   !mgf--
#else
   ! hydrophiliic BC
   real(r8), allocatable :: ss_alb_bc1     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_bc1    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_bc1(:) ! (numrad_snw);

   ! hydrophobic BC
   real(r8), allocatable :: ss_alb_bc2     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_bc2    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_bc2(:) ! (numrad_snw);
#endif

   ! hydrophobic OC
   real(r8), allocatable :: ss_alb_oc1     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_oc1    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_oc1(:) ! (numrad_snw);

   ! hydrophilic OC
   real(r8), allocatable :: ss_alb_oc2     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_oc2    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_oc2(:) ! (numrad_snw);

   ! dust species 1:
   real(r8), allocatable :: ss_alb_dst1     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_dst1    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_dst1(:) ! (numrad_snw);

   ! dust species 2:
   real(r8), allocatable :: ss_alb_dst2     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_dst2    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_dst2(:) ! (numrad_snw);

   ! dust species 3:
   real(r8), allocatable :: ss_alb_dst3     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_dst3    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_dst3(:) ! (numrad_snw);

   ! dust species 4:
   real(r8), allocatable :: ss_alb_dst4     (:) ! (numrad_snw);
   real(r8), allocatable :: asm_prm_dst4    (:) ! (numrad_snw);
   real(r8), allocatable :: ext_cff_mss_dst4(:) ! (numrad_snw);

#ifdef MODAL_AER
   ! Absorption enhancement factors for within-ice BC
   real(r8), allocatable :: bcenh (:,:,:) ! (numrad_snw,idx_bc_nclrds_max,idx_bcint_icerds_max);
#endif

   ! best-fit parameters for snow aging defined over:
   !  11 temperatures from 225 to 273 K
   !  31 temperature gradients from 0 to 300 K/m
   !   8 snow densities from 0 to 350 kg/m3
   ! (arrays declared here, but are set in iniTimeConst)
   !
   real(r8), allocatable :: snowage_tau   (:,:,:) ! (idx_rhos_max,idx_Tgrd_max,idx_T_max)  ! snow aging parameter retrieved from lookup table [hour]
   real(r8), allocatable :: snowage_kappa (:,:,:) ! (idx_rhos_max,idx_Tgrd_max,idx_T_max)  ! snow aging parameter retrieved from lookup table [unitless]
   real(r8), allocatable :: snowage_drdt0 (:,:,:) ! (idx_rhos_max,idx_Tgrd_max,idx_T_max)  ! snow aging parameter retrieved from lookup table [um hr-1]

   !
   ! !REVISION HISTORY:
   ! Created by Mark Flanner
   !-----------------------------------------------------------------------

CONTAINS

   SUBROUTINE SNICAR_RT (flg_snw_ice, flg_slr_in, &
                         coszen, snl, h2osno, frac_sno, &
                         h2osno_liq, h2osno_ice, snw_rds,   &
                         mss_cnc_aer_in, albsfc, albout, flx_abs)
!-----------------------------------------------------------------------
! !DESCRIPTION:
!  Determine reflectance of, and vertically-resolved solar absorption in,
!  snow with impurities.
!
!  Original references on physical models of snow reflectance include:
!  Wiscombe and Warren [1980] and Warren and Wiscombe [1980],
!  Journal of Atmospheric Sciences, 37,
!
!  The multi-layer solution for multiple-scattering used here is from:
!  Toon et al. [1989], Rapid calculation of radiative heating rates
!  and photodissociation rates in inhomogeneous multiple scattering atmospheres,
!  J. Geophys. Res., 94, D13, 16287-16301
!
!  The implementation of the SNICAR model in CLM/CSIM is described in:
!  Flanner, M., C. Zender, J. Randerson, and P. Rasch [2007],
!  Present-day climate forcing and response from black carbon in snow,
!  J. Geophys. Res., 112, D11202, doi: 10.1029/2006JD008003
!
! !USES:
!
!-----------------------------------------------------------------------
! !ARGUMENTS:

   IMPLICIT NONE

   integer  , intent(in)  :: flg_snw_ice                   ! flag: =1 when called from CLM, =2 when called from CSIM
   integer  , intent(in)  :: flg_slr_in                    ! flag: =1 for direct-beam incident flux,=2 for diffuse incident flux
   real(r8) , intent(in)  :: coszen                        ! cosine of solar zenith angle for next time step (col) [unitless]

   integer  , intent(in)  :: snl                           ! negative number of snow layers (col) [nbr]
   real(r8) , intent(in)  :: h2osno                        ! snow liquid water equivalent (col) [kg/m2]
   real(r8) , intent(in)  :: frac_sno                      ! fraction of ground covered by snow (0 to 1)

   real(r8) , intent(in)  :: h2osno_liq     ( maxsnl+1:0 ) ! liquid water content (col,lyr) [kg/m2]
   real(r8) , intent(in)  :: h2osno_ice     ( maxsnl+1:0 ) ! ice content (col,lyr) [kg/m2]
   integer  , intent(in)  :: snw_rds        ( maxsnl+1:0 ) ! snow effective radius (col,lyr) [microns, m^-6]
   real(r8) , intent(in)  :: mss_cnc_aer_in ( maxsnl+1:0 , 1:sno_nbr_aer ) ! mass concentration of all aerosol species (col,lyr,aer) [kg/kg]
   real(r8) , intent(in)  :: albsfc         ( 1:numrad )   ! albedo of surface underlying snow (col,bnd) [frc]
   real(r8) , intent(out) :: albout         ( 1:numrad )   ! snow albedo, averaged into 2 bands (=0 IF no sun or no snow) (col,bnd) [frc]
   real(r8) , intent(out) :: flx_abs        ( maxsnl+1:1 , 1:numrad ) ! absorbed flux in each layer per unit flux incident (col, lyr, bnd)
   !
   ! !LOCAL VARIABLES:
   !
   ! variables for snow radiative transfer calculations

   ! Local variables representing single-column values of arrays:
   integer :: snl_lcl                                   ! negative number of snow layers [nbr]
   integer :: snw_rds_lcl(maxsnl+1:0)                   ! snow effective radius [m^-6]
   real(r8):: flx_slrd_lcl(1:numrad_snw)                ! direct beam incident irradiance [W/m2] (set to 1)
   real(r8):: flx_slri_lcl(1:numrad_snw)                ! diffuse incident irradiance [W/m2] (set to 1)
   real(r8):: mss_cnc_aer_lcl(maxsnl+1:0,1:sno_nbr_aer) ! aerosol mass concentration (lyr,aer_nbr) [kg/kg]
   real(r8):: h2osno_lcl                                ! total column snow mass [kg/m2]
   real(r8):: h2osno_liq_lcl(maxsnl+1:0)                ! liquid water mass [kg/m2]
   real(r8):: h2osno_ice_lcl(maxsnl+1:0)                ! ice mass [kg/m2]
   real(r8):: albsfc_lcl(1:numrad_snw)                  ! albedo of underlying surface [frc]
   real(r8):: ss_alb_snw_lcl(maxsnl+1:0)                ! single-scatter albedo of ice grains (lyr) [frc]
   real(r8):: asm_prm_snw_lcl(maxsnl+1:0)               ! asymmetry parameter of ice grains (lyr) [frc]
   real(r8):: ext_cff_mss_snw_lcl(maxsnl+1:0)           ! mass extinction coefficient of ice grains (lyr) [m2/kg]
   real(r8):: ss_alb_aer_lcl(sno_nbr_aer)               ! single-scatter albedo of aerosol species (aer_nbr) [frc]
   real(r8):: asm_prm_aer_lcl(sno_nbr_aer)              ! asymmetry parameter of aerosol species (aer_nbr) [frc]
   real(r8):: ext_cff_mss_aer_lcl(sno_nbr_aer)          ! mass extinction coefficient of aerosol species (aer_nbr) [m2/kg]

#ifdef MODAL_AER
   !mgf++
   real(r8) :: rds_bcint_lcl(maxsnl+1:0)          ! effective radius of within-ice BC [nm]
   real(r8) :: rds_bcext_lcl(maxsnl+1:0)          ! effective radius of external BC [nm]
   !mgf--
#endif


   ! Other local variables
   integer :: APRX_TYP                            ! two-stream approximation type
                                                  ! (1=Eddington, 2=Quadrature, 3=Hemispheric Mean) [nbr]
   integer :: DELTA                               ! flag to USE Delta approximation (Joseph, 1976)
                                                  ! (1= USE, 0= don't USE)
   real(r8):: flx_wgt(1:numrad_snw)               ! weights applied to spectral bands,
                                                  ! specific to direct and diffuse cases (bnd) [frc]

   integer :: flg_nosnl                           ! flag: =1 IF there is snow, but zero snow layers,
                                                  ! =0 IF at least 1 snow layer [flg]
   integer :: trip                                ! flag: =1 to redo RT calculation IF result is unrealistic
   integer :: flg_dover                           ! defines conditions for RT redo (explained below)

   real(r8):: albedo                              ! temporary snow albedo [frc]
   real(r8):: flx_sum                             ! temporary summation variable for NIR weighting
   real(r8):: albout_lcl(numrad_snw)              ! snow albedo by band [frc]
   real(r8):: flx_abs_lcl(maxsnl+1:1,numrad_snw)  ! absorbed flux per unit incident flux at top of snowpack (lyr,bnd) [frc]

   real(r8):: L_snw(maxsnl+1:0)                   ! h2o mass (liquid+solid) in snow layer (lyr) [kg/m2]
   real(r8):: tau_snw(maxsnl+1:0)                 ! snow optical depth (lyr) [unitless]
   real(r8):: L_aer(maxsnl+1:0,sno_nbr_aer)       ! aerosol mass in snow layer (lyr,nbr_aer) [kg/m2]
   real(r8):: tau_aer(maxsnl+1:0,sno_nbr_aer)     ! aerosol optical depth (lyr,nbr_aer) [unitless]
   real(r8):: tau_sum                             ! cumulative (snow+aerosol) optical depth [unitless]
   real(r8):: tau_elm(maxsnl+1:0)                 ! column optical depth from layer bottom to snowpack top (lyr) [unitless]
   real(r8):: omega_sum                           ! temporary summation of single-scatter albedo of all aerosols [frc]
   real(r8):: g_sum                               ! temporary summation of asymmetry parameter of all aerosols [frc]

   real(r8):: tau(maxsnl+1:0)                     ! weighted optical depth of snow+aerosol layer (lyr) [unitless]
   real(r8):: omega(maxsnl+1:0)                   ! weighted single-scatter albedo of snow+aerosol layer (lyr) [frc]
   real(r8):: g(maxsnl+1:0)                       ! weighted asymmetry parameter of snow+aerosol layer (lyr) [frc]
   real(r8):: tau_star(maxsnl+1:0)                ! transformed (i.e. Delta-Eddington) optical depth of snow+aerosol layer
                                                  ! (lyr) [unitless]
   real(r8):: omega_star(maxsnl+1:0)              ! transformed (i.e. Delta-Eddington) SSA of snow+aerosol layer (lyr) [frc]
   real(r8):: g_star(maxsnl+1:0)                  ! transformed (i.e. Delta-Eddington) asymmetry paramater of snow+aerosol layer
                                                  ! (lyr) [frc]

   integer :: bnd_idx                             ! spectral band index (1 <= bnd_idx <= numrad_snw) [idx]
   integer :: rds_idx                             ! snow effective radius index for retrieving
                                                  ! Mie parameters from lookup table [idx]
   integer :: snl_btm                             ! index of bottom snow layer (0) [idx]
   integer :: snl_top                             ! index of top snow layer (-4 to 0) [idx]
   integer :: fc                                  ! column filter index
   integer :: i                                   ! layer index [idx]
   integer :: j                                   ! aerosol number index [idx]
   integer :: n                                   ! tridiagonal matrix index [idx]
   integer :: m                                   ! secondary layer index [idx]
   integer :: nint_snw_rds_min                    ! nearest integer value of snw_rds_min

   real(r8):: F_direct(maxsnl+1:0)                ! direct-beam radiation at bottom of layer interface (lyr) [W/m^2]
   real(r8):: F_net(maxsnl+1:0)                   ! net radiative flux at bottom of layer interface (lyr) [W/m^2]
   real(r8):: F_abs(maxsnl+1:0)                   ! net absorbed radiative energy (lyr) [W/m^2]
   real(r8):: F_abs_sum                           ! total absorbed energy in column [W/m^2]
   real(r8):: F_sfc_pls                           ! upward radiative flux at snowpack top [W/m^2]
   real(r8):: F_btm_net                           ! net flux at bottom of snowpack [W/m^2]
   real(r8):: F_sfc_net                           ! net flux at top of snowpack [W/m^2]
   real(r8):: energy_sum                          ! sum of all energy terms; should be 0.0 [W/m^2]
   real(r8):: F_direct_btm                        ! direct-beam radiation at bottom of snowpack [W/m^2]
   real(r8):: mu_not                              ! cosine of solar zenith angle (used locally) [frc]

   integer :: err_idx                             ! counter for number of times through error loop [nbr]
   real(r8):: pi                                  ! 3.1415...

   ! intermediate variables for radiative transfer approximation:
   real(r8):: gamma1(maxsnl+1:0)    ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: gamma2(maxsnl+1:0)    ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: gamma3(maxsnl+1:0)    ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: gamma4(maxsnl+1:0)    ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: lambda(maxsnl+1:0)    ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: GAMMA(maxsnl+1:0)     ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: mu_one                ! two-stream coefficient from Toon et al. (lyr) [unitless]
   real(r8):: e1(maxsnl+1:0)        ! tri-diag intermediate variable from Toon et al. (lyr)
   real(r8):: e2(maxsnl+1:0)        ! tri-diag intermediate variable from Toon et al. (lyr)
   real(r8):: e3(maxsnl+1:0)        ! tri-diag intermediate variable from Toon et al. (lyr)
   real(r8):: e4(maxsnl+1:0)        ! tri-diag intermediate variable from Toon et al. (lyr)
   real(r8):: C_pls_btm(maxsnl+1:0) ! intermediate variable: upward flux at bottom interface (lyr) [W/m2]
   real(r8):: C_mns_btm(maxsnl+1:0) ! intermediate variable: downward flux at bottom interface (lyr) [W/m2]
   real(r8):: C_pls_top(maxsnl+1:0) ! intermediate variable: upward flux at top interface (lyr) [W/m2]
   real(r8):: C_mns_top(maxsnl+1:0) ! intermediate variable: downward flux at top interface (lyr) [W/m2]
   real(r8):: A(2*maxsnl+1:0)       ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: B(2*maxsnl+1:0)       ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: D(2*maxsnl+1:0)       ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: E(2*maxsnl+1:0)       ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: AS(2*maxsnl+1:0)      ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: DS(2*maxsnl+1:0)      ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: X(2*maxsnl+1:0)       ! tri-diag intermediate variable from Toon et al. (2*lyr)
   real(r8):: Y(2*maxsnl+1:0)       ! tri-diag intermediate variable from Toon et al. (2*lyr)
   !-----------------------------------------------------------------------
#ifdef MODAL_AER
   !mgf++
   integer :: idx_bcint_icerds                   ! index of ice effective radius for optical properties lookup table
   integer :: idx_bcint_nclrds                   ! index of within-ice BC effective radius for optical properties lookup table
   integer :: idx_bcext_nclrds                   ! index of external BC effective radius for optical properties lookup table
   real(r8):: enh_fct                            ! extinction/absorption enhancement factor for within-ice BC
   real(r8):: tmp1                               ! temporary variable
   !mgf--
#endif

   ! Enforce expected array sizes

 !  associate(&
 !       snl      => col_pp%snl               , & ! Input:  [integer (:)]  negative number of snow layers (col) [nbr]
 !       h2osno   => col_ws%h2osno            , & ! Input:  [real(r8) (:)]  snow liquid water equivalent (col) [kg/m2]
 !       frac_sno =>  col_ws%frac_sno_eff       & ! Input:  [real(r8) (:)]  fraction of ground covered by snow (0 to 1)
 !       )

      ! Define constants
      pi = SHR_CONST_PI
      nint_snw_rds_min = nint(snw_rds_min)

      ! always USE Delta approximation for snow
      DELTA = 1

      ! (when called from CSIM, there is only one column)

      ! Zero absorbed radiative fluxes:
      DO i=maxsnl+1,1,1
         flx_abs_lcl(:,:)   = 0._r8
         flx_abs(i,:) = 0._r8
      ENDDO

      ! set snow/ice mass to be used for RT:
      IF (flg_snw_ice == 1) THEN
         h2osno_lcl = h2osno
      ELSE
         h2osno_lcl = h2osno_ice(0)
      ENDIF

      ! Qualifier for computing snow RT:
      !  1) sunlight from atmosphere model
      !  2) minimum amount of snow on ground.
      !     Otherwise, set snow albedo to zero
      IF ((coszen > 0._r8) .and. (h2osno_lcl > min_snw)) THEN

         ! Set variables specific to CLM
         IF (flg_snw_ice == 1) THEN
            ! If there is snow, but zero snow layers, we must create a layer locally.
            ! This layer is presumed to have the fresh snow effective radius.
            IF (snl > -1) THEN
               flg_nosnl         =  1
               snl_lcl           =  -1
               h2osno_ice_lcl(0) =  h2osno_lcl
               h2osno_liq_lcl(0) =  0._r8
               snw_rds_lcl(0)    =  nint_snw_rds_min
            ELSE
               flg_nosnl         =  0
               snl_lcl           =  snl
               h2osno_liq_lcl(:) =  h2osno_liq(:)
               h2osno_ice_lcl(:) =  h2osno_ice(:)
               snw_rds_lcl(:)    =  snw_rds(:)
            ENDIF

            snl_btm   = 0
            snl_top   = snl_lcl+1

            ! Set variables specific to CSIM
         ELSE
            flg_nosnl         = 0
            snl_lcl           = -1
            h2osno_liq_lcl(:) = h2osno_liq(:)
            h2osno_ice_lcl(:) = h2osno_ice(:)
            snw_rds_lcl(:)    = snw_rds(:)
            snl_btm           = 0
            snl_top           = 0
         ENDIF

#ifdef MODAL_AER
         !mgf++
         !
         ! Assume fixed BC effective radii of 100nm. This is close to
         ! the effective radius of 95nm (number median radius of
         ! 40nm) assumed for freshly-emitted BC in MAM.  Future
         ! implementations may prognose the BC effective radius in
         ! snow.
         rds_bcint_lcl(:)  =  100._r8
         rds_bcext_lcl(:)  =  100._r8
         !mgf--
#endif

         ! Set local aerosol array
         DO j=1,sno_nbr_aer
            mss_cnc_aer_lcl(:,j) = mss_cnc_aer_in(:,j)
         ENDDO


         ! Set spectral underlying surface albedos to their corresponding VIS or NIR albedos
         albsfc_lcl(1)                       = albsfc(1)
         albsfc_lcl(nir_bnd_bgn:nir_bnd_end) = albsfc(2)


         ! Error check for snow grain size:
#ifndef _OPENACC
         IF (p_is_master) THEN
            DO i=snl_top,snl_btm,1
               IF ((snw_rds_lcl(i) < snw_rds_min_tbl) .or. (snw_rds_lcl(i) > snw_rds_max_tbl)) THEN
                  write (iulog,*)  "SNICAR ERROR: snow grain radius of out of bounds."
                  write (iulog,*) "flg_snw_ice= ", flg_snw_ice
                  write (iulog,*) " level: ", i, " snl(c)= ", snl_lcl
                  write (iulog,*) "h2osno(c)= ", h2osno_lcl
                  CALL abort
               ENDIF
            ENDDO
         ENDIF
#endif

         ! Incident flux weighting parameters
         !  - sum of all VIS bands must equal 1
         !  - sum of all NIR bands must equal 1
         !
         ! Spectral bands (5-band CASE)
         !  Band 1: 0.3-0.7um (VIS)
         !  Band 2: 0.7-1.0um (NIR)
         !  Band 3: 1.0-1.2um (NIR)
         !  Band 4: 1.2-1.5um (NIR)
         !  Band 5: 1.5-5.0um (NIR)
         !
         ! The following weights are appropriate for surface-incident flux in a mid-latitude winter atmosphere
         !
         ! 3-band weights
         IF (numrad_snw==3) THEN
            ! Direct:
            IF (flg_slr_in == 1) THEN
               flx_wgt(1) = 1._r8
               flx_wgt(2) = 0.66628670195247_r8
               flx_wgt(3) = 0.33371329804753_r8
               ! Diffuse:
            ELSEIF (flg_slr_in == 2) THEN
               flx_wgt(1) = 1._r8
               flx_wgt(2) = 0.77887652162877_r8
               flx_wgt(3) = 0.22112347837123_r8
            ENDIF

            ! 5-band weights
         ELSEIF(numrad_snw==5) THEN
            ! Direct:
            IF (flg_slr_in == 1) THEN
               flx_wgt(1) = 1._r8
               flx_wgt(2) = 0.49352158521175_r8
               flx_wgt(3) = 0.18099494230665_r8
               flx_wgt(4) = 0.12094898498813_r8
               flx_wgt(5) = 0.20453448749347_r8
               ! Diffuse:
            ELSEIF (flg_slr_in == 2) THEN
               flx_wgt(1) = 1._r8
               flx_wgt(2) = 0.58581507618433_r8
               flx_wgt(3) = 0.20156903770812_r8
               flx_wgt(4) = 0.10917889346386_r8
               flx_wgt(5) = 0.10343699264369_r8
            ENDIF
         ENDIF

         ! Loop over snow spectral bands
         DO bnd_idx = 1,numrad_snw

            mu_not    = coszen  ! must set here, because of error handling
            flg_dover = 1              ! default is to redo
            err_idx   = 0              ! number of times through loop

            DO WHILE (flg_dover > 0)

               ! DEFAULT APPROXIMATIONS:
               !  VIS:       Delta-Eddington
               !  NIR (all): Delta-Hemispheric Mean
               !  WARNING:   DO NOT USE DELTA-EDDINGTON FOR NIR DIFFUSE - this sometimes results in negative albedo
               !
               ! ERROR CONDITIONS:
               !  Conditions which cause "trip", resulting in redo of RT approximation:
               !   1. negative absorbed flux
               !   2. total absorbed flux greater than incident flux
               !   3. negative albedo
               !   NOTE: These errors have only been encountered in spectral bands 4 and 5
               !
               ! ERROR HANDLING
               !  1st error (flg_dover=2): switch approximation (Edd->HM or HM->Edd)
               !  2nd error (flg_dover=3): change zenith angle by 0.02 (this happens about 1 in 10^6 cases)
               !  3rd error (flg_dover=4): switch approximation with new zenith
               !  Subsequent errors: repeatedly change zenith and approximations...

               IF (bnd_idx == 1) THEN
                  IF (flg_dover == 2) THEN
                     APRX_TYP = 3
                  ELSEIF (flg_dover == 3) THEN
                     APRX_TYP = 1
                     IF (coszen > 0.5_r8) THEN
                        mu_not = mu_not - 0.02_r8
                     ELSE
                        mu_not = mu_not + 0.02_r8
                     ENDIF
                  ELSEIF (flg_dover == 4) THEN
                     APRX_TYP = 3
                  ELSE
                     APRX_TYP = 1
                  ENDIF

               ELSE
                  IF (flg_dover == 2) THEN
                     APRX_TYP = 1
                  ELSEIF (flg_dover == 3) THEN
                     APRX_TYP = 3
                     IF (coszen > 0.5_r8) THEN
                        mu_not = mu_not - 0.02_r8
                     ELSE
                        mu_not = mu_not + 0.02_r8
                     ENDIF
                  ELSEIF (flg_dover == 4) THEN
                     APRX_TYP = 1
                  ELSE
                     APRX_TYP = 3
                  ENDIF

               ENDIF

               ! Set direct or diffuse incident irradiance to 1
               ! (This has to be within the bnd loop because mu_not is adjusted in rare cases)
               IF (flg_slr_in == 1) THEN
                  flx_slrd_lcl(bnd_idx) = 1._r8/(mu_not*pi) ! this corresponds to incident irradiance of 1.0
                  flx_slri_lcl(bnd_idx) = 0._r8
               ELSE
                  flx_slrd_lcl(bnd_idx) = 0._r8
                  flx_slri_lcl(bnd_idx) = 1._r8
               ENDIF

               ! Pre-emptive error handling: aerosols can reap havoc on these absorptive bands.
               ! Since extremely high soot concentrations have a negligible effect on these bands, zero them.
               IF ( (numrad_snw == 5).and.((bnd_idx == 5).or.(bnd_idx == 4)) ) THEN
                  mss_cnc_aer_lcl(:,:) = 0._r8
               ENDIF

               IF ( (numrad_snw == 3).and.(bnd_idx == 3) ) THEN
                  mss_cnc_aer_lcl(:,:) = 0._r8
               ENDIF

               ! Define local Mie parameters based on snow grain size and aerosol species,
               !  retrieved from a lookup table.
               IF (flg_slr_in == 1) THEN
                  DO i=snl_top,snl_btm,1
                     rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
                     ! snow optical properties (direct radiation)
                     ss_alb_snw_lcl(i)      = ss_alb_snw_drc(rds_idx,bnd_idx)
                     asm_prm_snw_lcl(i)     = asm_prm_snw_drc(rds_idx,bnd_idx)
                     ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_drc(rds_idx,bnd_idx)
                  ENDDO
               ELSEIF (flg_slr_in == 2) THEN
                  DO i=snl_top,snl_btm,1
                     rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
                     ! snow optical properties (diffuse radiation)
                     ss_alb_snw_lcl(i)      = ss_alb_snw_dfs(rds_idx,bnd_idx)
                     asm_prm_snw_lcl(i)     = asm_prm_snw_dfs(rds_idx,bnd_idx)
                     ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_dfs(rds_idx,bnd_idx)
                  ENDDO
               ENDIF

!H. Wang
              ! aerosol species 1 optical properties
              ! ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx)
              ! asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx)
              ! ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx)

              ! aerosol species 2 optical properties
              ! ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx)
              ! asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx)
              ! ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx)
!H. Wang
              ! aerosol species 3 optical properties
               ss_alb_aer_lcl(3)        = ss_alb_oc1(bnd_idx)
               asm_prm_aer_lcl(3)       = asm_prm_oc1(bnd_idx)
               ext_cff_mss_aer_lcl(3)   = ext_cff_mss_oc1(bnd_idx)

              ! aerosol species 4 optical properties
               ss_alb_aer_lcl(4)        = ss_alb_oc2(bnd_idx)
               asm_prm_aer_lcl(4)       = asm_prm_oc2(bnd_idx)
               ext_cff_mss_aer_lcl(4)   = ext_cff_mss_oc2(bnd_idx)

              ! aerosol species 5 optical properties
               ss_alb_aer_lcl(5)        = ss_alb_dst1(bnd_idx)
               asm_prm_aer_lcl(5)       = asm_prm_dst1(bnd_idx)
               ext_cff_mss_aer_lcl(5)   = ext_cff_mss_dst1(bnd_idx)

              ! aerosol species 6 optical properties
               ss_alb_aer_lcl(6)        = ss_alb_dst2(bnd_idx)
               asm_prm_aer_lcl(6)       = asm_prm_dst2(bnd_idx)
               ext_cff_mss_aer_lcl(6)   = ext_cff_mss_dst2(bnd_idx)

              ! aerosol species 7 optical properties
               ss_alb_aer_lcl(7)        = ss_alb_dst3(bnd_idx)
               asm_prm_aer_lcl(7)       = asm_prm_dst3(bnd_idx)
               ext_cff_mss_aer_lcl(7)   = ext_cff_mss_dst3(bnd_idx)

              ! aerosol species 8 optical properties
               ss_alb_aer_lcl(8)        = ss_alb_dst4(bnd_idx)
               asm_prm_aer_lcl(8)       = asm_prm_dst4(bnd_idx)
               ext_cff_mss_aer_lcl(8)   = ext_cff_mss_dst4(bnd_idx)


               ! 1. snow and aerosol layer column mass (L_snw, L_aer [kg/m^2])
               ! 2. optical Depths (tau_snw, tau_aer)
               ! 3. weighted Mie properties (tau, omega, g)

               ! Weighted Mie parameters of each layer
               DO i=snl_top,snl_btm,1
#ifdef MODAL_AER
                 !mgf++ within-ice and external BC optical properties
                 !
                 ! Lookup table indices for BC optical properties,
                 ! dependent on snow grain size and BC particle
                 ! size.

                 ! valid for 25 < snw_rds < 1625 um:
                 IF (snw_rds_lcl(i) < 125) THEN
                    tmp1 = snw_rds_lcl(i)/50
                    idx_bcint_icerds = nint(tmp1)
                 ELSEIF (snw_rds_lcl(i) < 175) THEN
                    idx_bcint_icerds = 2
                 ELSE
                    tmp1 = (snw_rds_lcl(i)/250)+2
                    idx_bcint_icerds = nint(tmp1)
                 ENDIF

                 ! valid for 25 < bc_rds < 525 nm
                 idx_bcint_nclrds = nint(rds_bcint_lcl(i)/50)
                 idx_bcext_nclrds = nint(rds_bcext_lcl(i)/50)

                 ! check bounds:
                 IF (idx_bcint_icerds < idx_bcint_icerds_min) idx_bcint_icerds = idx_bcint_icerds_min
                 IF (idx_bcint_icerds > idx_bcint_icerds_max) idx_bcint_icerds = idx_bcint_icerds_max
                 IF (idx_bcint_nclrds < idx_bc_nclrds_min) idx_bcint_nclrds = idx_bc_nclrds_min
                 IF (idx_bcint_nclrds > idx_bc_nclrds_max) idx_bcint_nclrds = idx_bc_nclrds_max
                 IF (idx_bcext_nclrds < idx_bc_nclrds_min) idx_bcext_nclrds = idx_bc_nclrds_min
                 IF (idx_bcext_nclrds > idx_bc_nclrds_max) idx_bcext_nclrds = idx_bc_nclrds_max

                 ! retrieve absorption enhancement factor for within-ice BC
                 enh_fct = bcenh(bnd_idx,idx_bcint_nclrds,idx_bcint_icerds)

                 ! get BC optical properties (moved from above)
                 ! aerosol species 1 optical properties (within-ice BC)
                 ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx,idx_bcint_nclrds)
                 asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx,idx_bcint_nclrds)
                 ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx,idx_bcint_nclrds)*enh_fct

                 ! aerosol species 2 optical properties (external BC)
                 ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx,idx_bcext_nclrds)
                 asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx,idx_bcext_nclrds)
                 ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx,idx_bcext_nclrds)

#else
                 ! bulk aerosol treatment (BC optical properties independent
                 ! of BC and ice grain size)
                 ! aerosol species 1 optical properties (within-ice BC)
                 ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx)
                 asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx)
                 ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx)

                 ! aerosol species 2 optical properties
                 ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx)
                 asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx)
                 ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx)
#endif
                 !mgf--

                 L_snw(i)   = h2osno_ice_lcl(i)+h2osno_liq_lcl(i)
                 tau_snw(i) = L_snw(i)*ext_cff_mss_snw_lcl(i)

                 DO j=1,sno_nbr_aer
                    L_aer(i,j)   = L_snw(i)*mss_cnc_aer_lcl(i,j)
                    tau_aer(i,j) = L_aer(i,j)*ext_cff_mss_aer_lcl(j)
                 ENDDO

                 tau_sum   = 0._r8
                 omega_sum = 0._r8
                 g_sum     = 0._r8

                 DO j=1,sno_nbr_aer
                    tau_sum    = tau_sum + tau_aer(i,j)
                    omega_sum  = omega_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j))
                    g_sum      = g_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j)*asm_prm_aer_lcl(j))
                 ENDDO

                 tau(i)    = tau_sum + tau_snw(i)
                 omega(i)  = (1/tau(i))*(omega_sum+(ss_alb_snw_lcl(i)*tau_snw(i)))
                 g(i)      = (1/(tau(i)*omega(i)))*(g_sum+ (asm_prm_snw_lcl(i)*ss_alb_snw_lcl(i)*tau_snw(i)))
               ENDDO

               ! DELTA transformations, IF requested
               IF (DELTA == 1) THEN
                  DO i=snl_top,snl_btm,1
                     g_star(i)     = g(i)/(1+g(i))
                     omega_star(i) = ((1-(g(i)**2))*omega(i)) / (1-(omega(i)*(g(i)**2)))
                     tau_star(i)   = (1-(omega(i)*(g(i)**2)))*tau(i)
                  ENDDO
               ELSE
                  DO i=snl_top,snl_btm,1
                     g_star(i)     = g(i)
                     omega_star(i) = omega(i)
                     tau_star(i)   = tau(i)
                  ENDDO
               ENDIF

               ! Total column optical depth:
               ! tau_elm(i) = total optical depth above the bottom of layer i
               tau_elm(snl_top) = 0._r8
               DO i=snl_top+1,snl_btm,1
                  tau_elm(i) = tau_elm(i-1)+tau_star(i-1)
               ENDDO

               ! Direct radiation at bottom of snowpack:
               F_direct_btm = albsfc_lcl(bnd_idx)*mu_not * &
                    exp(-(tau_elm(snl_btm)+tau_star(snl_btm))/mu_not)*pi*flx_slrd_lcl(bnd_idx)

               ! Intermediates
               ! Gamma values are approximation-specific.

               ! Eddington
               IF (APRX_TYP==1) THEN
                  DO i=snl_top,snl_btm,1
                     gamma1(i) = (7-(omega_star(i)*(4+(3*g_star(i)))))/4
                     gamma2(i) = -(1-(omega_star(i)*(4-(3*g_star(i)))))/4
                     gamma3(i) = (2-(3*g_star(i)*mu_not))/4
                     gamma4(i) = 1-gamma3(i)
                     mu_one    = 0.5
                  ENDDO

               ! Quadrature
               ELSEIF (APRX_TYP==2) THEN
                  DO i=snl_top,snl_btm,1
                     gamma1(i) = (3**0.5)*(2-(omega_star(i)*(1+g_star(i))))/2
                     gamma2(i) = omega_star(i)*(3**0.5)*(1-g_star(i))/2
                     gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2
                     gamma4(i) = 1-gamma3(i)
                     mu_one    = 1/(3**0.5)
                  ENDDO

               ! Hemispheric Mean
               ELSEIF (APRX_TYP==3) THEN
                  DO i=snl_top,snl_btm,1
                     gamma1(i) = 2 - (omega_star(i)*(1+g_star(i)))
                     gamma2(i) = omega_star(i)*(1-g_star(i))
                     gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2
                     gamma4(i) = 1-gamma3(i)
                     mu_one    = 0.5
                  ENDDO
               ENDIF

               ! Intermediates for tri-diagonal solution
               DO i=snl_top,snl_btm,1
                  lambda(i) = sqrt(abs((gamma1(i)**2) - (gamma2(i)**2)))
                  GAMMA(i)  = gamma2(i)/(gamma1(i)+lambda(i))

                  e1(i)     = 1+(GAMMA(i)*exp(-lambda(i)*tau_star(i)))
                  e2(i)     = 1-(GAMMA(i)*exp(-lambda(i)*tau_star(i)))
                  e3(i)     = GAMMA(i) + exp(-lambda(i)*tau_star(i))
                  e4(i)     = GAMMA(i) - exp(-lambda(i)*tau_star(i))
               ENDDO !ENDDO over snow layers

               ! Intermediates for tri-diagonal solution
               DO i=snl_top,snl_btm,1
                  IF (flg_slr_in == 1) THEN

                     C_pls_btm(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
                          exp(-(tau_elm(i)+tau_star(i))/mu_not)*   &
                          (((gamma1(i)-(1/mu_not))*gamma3(i))+     &
                          (gamma4(i)*gamma2(i))))/((lambda(i)**2)-(1/(mu_not**2)))

                     C_mns_btm(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
                          exp(-(tau_elm(i)+tau_star(i))/mu_not)*   &
                          (((gamma1(i)+(1/mu_not))*gamma4(i))+     &
                          (gamma2(i)*gamma3(i))))/((lambda(i)**2)-(1/(mu_not**2)))

                     C_pls_top(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
                          exp(-tau_elm(i)/mu_not)*(((gamma1(i)-(1/mu_not))* &
                          gamma3(i))+(gamma4(i)*gamma2(i))))/((lambda(i)**2)-(1/(mu_not**2)))

                     C_mns_top(i) = (omega_star(i)*pi*flx_slrd_lcl(bnd_idx)* &
                          exp(-tau_elm(i)/mu_not)*(((gamma1(i)+(1/mu_not))* &
                          gamma4(i))+(gamma2(i)*gamma3(i))))/((lambda(i)**2)-(1/(mu_not**2)))

                  ELSE
                     C_pls_btm(i) = 0._r8
                     C_mns_btm(i) = 0._r8
                     C_pls_top(i) = 0._r8
                     C_mns_top(i) = 0._r8
                  ENDIF
               ENDDO

               ! Coefficients for tridiaganol matrix solution
               DO i=2*snl_lcl+1,0,1

                  !Boundary values for i=1 and i=2*snl_lcl, specifics for i=odd and i=even
                  IF (i==(2*snl_lcl+1)) THEN
                     A(i) = 0
                     B(i) = e1(snl_top)
                     D(i) = -e2(snl_top)
                     E(i) = flx_slri_lcl(bnd_idx)-C_mns_top(snl_top)

                  ELSEIF(i==0) THEN
                     A(i) = e1(snl_btm)-(albsfc_lcl(bnd_idx)*e3(snl_btm))
                     B(i) = e2(snl_btm)-(albsfc_lcl(bnd_idx)*e4(snl_btm))
                     D(i) = 0
                     E(i) = F_direct_btm-C_pls_btm(snl_btm)+(albsfc_lcl(bnd_idx)*C_mns_btm(snl_btm))

                  ELSEIF(mod(i,2)==-1) THEN   ! If odd and i>=3 (n=1 for i=3)
                     n=floor(i/2.0)
                     A(i) = (e2(n)*e3(n))-(e4(n)*e1(n))
                     B(i) = (e1(n)*e1(n+1))-(e3(n)*e3(n+1))
                     D(i) = (e3(n)*e4(n+1))-(e1(n)*e2(n+1))
                     E(i) = (e3(n)*(C_pls_top(n+1)-C_pls_btm(n)))+(e1(n)*(C_mns_btm(n)-C_mns_top(n+1)))

                  ELSEIF(mod(i,2)==0) THEN    ! If even and i<=2*snl_lcl
                     n=(i/2)
                     A(i) = (e2(n+1)*e1(n))-(e3(n)*e4(n+1))
                     B(i) = (e2(n)*e2(n+1))-(e4(n)*e4(n+1))
                     D(i) = (e1(n+1)*e4(n+1))-(e2(n+1)*e3(n+1))
                     E(i) = (e2(n+1)*(C_pls_top(n+1)-C_pls_btm(n)))+(e4(n+1)*(C_mns_top(n+1)-C_mns_btm(n)))
                  ENDIF
               ENDDO

               AS(0) = A(0)/B(0)
               DS(0) = E(0)/B(0)

               DO i=-1,(2*snl_lcl+1),-1
                  X(i)  = 1/(B(i)-(D(i)*AS(i+1)))
                  AS(i) = A(i)*X(i)
                  DS(i) = (E(i)-(D(i)*DS(i+1)))*X(i)
               ENDDO

               Y(2*snl_lcl+1) = DS(2*snl_lcl+1)
               DO i=(2*snl_lcl+2),0,1
                  Y(i) = DS(i)-(AS(i)*Y(i-1))
               ENDDO

               ! Downward direct-beam and net flux (F_net) at the base of each layer:
               DO i=snl_top,snl_btm,1
                  F_direct(i) = mu_not*pi*flx_slrd_lcl(bnd_idx)*exp(-(tau_elm(i)+tau_star(i))/mu_not)
                  F_net(i)    = (Y(2*i-1)*(e1(i)-e3(i))) + (Y(2*i)*(e2(i)-e4(i))) + &
                       C_pls_btm(i) - C_mns_btm(i) - F_direct(i)
               ENDDO

               ! Upward flux at snowpack top:
               F_sfc_pls = (Y(2*snl_lcl+1)*(exp(-lambda(snl_top)*tau_star(snl_top))+ &
                    GAMMA(snl_top))) + (Y(2*snl_lcl+2)*(exp(-lambda(snl_top)* &
                    tau_star(snl_top))-GAMMA(snl_top))) + C_pls_top(snl_top)

               ! Net flux at bottom = absorbed radiation by underlying surface:
               F_btm_net = -F_net(snl_btm)


               ! Bulk column albedo and surface net flux
               albedo    = F_sfc_pls/((mu_not*pi*flx_slrd_lcl(bnd_idx))+flx_slri_lcl(bnd_idx))
               F_sfc_net = F_sfc_pls - ((mu_not*pi*flx_slrd_lcl(bnd_idx))+flx_slri_lcl(bnd_idx))

               trip = 0
               ! Absorbed flux in each layer
               DO i=snl_top,snl_btm,1
                  IF(i==snl_top) THEN
                     F_abs(i) = F_net(i)-F_sfc_net
                  ELSE
                     F_abs(i) = F_net(i)-F_net(i-1)
                  ENDIF
                  flx_abs_lcl(i,bnd_idx) = F_abs(i)


                  ! ERROR check: negative absorption
                  IF (flx_abs_lcl(i,bnd_idx) < -0.00001) THEN
                     trip = 1
                  ENDIF
               ENDDO

               flx_abs_lcl(1,bnd_idx) = F_btm_net

               IF (flg_nosnl == 1) THEN
                  ! If there are no snow layers (but still snow), all absorbed energy must be in top soil layer
                  !flx_abs_lcl(:,bnd_idx) = 0._r8
                  !flx_abs_lcl(1,bnd_idx) = F_abs(0) + F_btm_net

                  ! changed on 20070408:
                  ! OK to put absorbed energy in the fictitous snow layer because routine SurfaceRadiation
                  ! handles the CASE of no snow layers. Then, IF a snow layer is addded between now and
                  ! SurfaceRadiation (called in CanopyHydrology), absorbed energy will be properly distributed.
                  flx_abs_lcl(0,bnd_idx) = F_abs(0)
                  flx_abs_lcl(1,bnd_idx) = F_btm_net

               ENDIF

               !Underflow check (we've already tripped the error condition above)
               DO i=snl_top,1,1
                  IF (flx_abs_lcl(i,bnd_idx) < 0._r8) THEN
                     flx_abs_lcl(i,bnd_idx) = 0._r8
                  ENDIF
               ENDDO

               F_abs_sum = 0._r8
               DO i=snl_top,snl_btm,1
                  F_abs_sum = F_abs_sum + F_abs(i)
               ENDDO


               !ERROR check: absorption greater than incident flux
               ! (should make condition more generic than "1._r8")
               IF (F_abs_sum > 1._r8) THEN
                  trip = 1
               ENDIF

               !ERROR check:
               IF ((albedo < 0._r8).and.(trip==0)) THEN
                  trip = 1
               ENDIF

               ! Set conditions for redoing RT calculation
               IF ((trip == 1).and.(flg_dover == 1)) THEN
                  flg_dover = 2
               ELSEIF ((trip == 1).and.(flg_dover == 2)) THEN
                  flg_dover = 3
               ELSEIF ((trip == 1).and.(flg_dover == 3)) THEN
                  flg_dover = 4
               ELSEIF((trip == 1).and.(flg_dover == 4).and.(err_idx < 20)) THEN
                  flg_dover = 3
                  err_idx = err_idx + 1
               ELSEIF((trip == 1).and.(flg_dover == 4).and.(err_idx >= 20)) THEN
                 flg_dover = 0
#ifndef _OPENACC
                 IF (p_is_master) THEN
                    write(iulog,*) "SNICAR ERROR: FOUND A WORMHOLE. STUCK IN INFINITE LOOP! Called from: ", flg_snw_ice
                    write(iulog,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(0)
                    write(iulog,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
                    write(iulog,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
                    write(iulog,*) "SNICAR STATS: soot1(0)= ", mss_cnc_aer_lcl(0,1)
                    write(iulog,*) "SNICAR STATS: soot2(0)= ", mss_cnc_aer_lcl(0,2)
                    write(iulog,*) "SNICAR STATS: dust1(0)= ", mss_cnc_aer_lcl(0,3)
                    write(iulog,*) "SNICAR STATS: dust2(0)= ", mss_cnc_aer_lcl(0,4)
                    write(iulog,*) "SNICAR STATS: dust3(0)= ", mss_cnc_aer_lcl(0,5)
                    write(iulog,*) "SNICAR STATS: dust4(0)= ", mss_cnc_aer_lcl(0,6)
                    write(iulog,*) "frac_sno: ", frac_sno
                    CALL abort
                 ENDIF
#endif
               ELSE
                  flg_dover = 0
               ENDIF

            ENDDO !ENDDO WHILE (flg_dover > 0)

            ! Energy conservation check:
            ! Incident direct+diffuse radiation equals (absorbed+bulk_transmitted+bulk_reflected)
            energy_sum = (mu_not*pi*flx_slrd_lcl(bnd_idx)) + flx_slri_lcl(bnd_idx) - (F_abs_sum + F_btm_net + F_sfc_pls)
            IF (abs(energy_sum) > 0.00001_r8) THEN
#ifndef _OPENACC
               IF (p_is_master) THEN
                  write(iulog,*) "SNICAR ERROR: Energy conservation error of : ", energy_sum
                  CALL abort
               ENDIF
#endif
            ENDIF

            albout_lcl(bnd_idx) = albedo

            ! Check that albedo is less than 1
            IF (albout_lcl(bnd_idx) > 1.0) THEN
#ifndef _OPENACC
               IF (p_is_master) THEN
                  write(iulog,*) "SNICAR ERROR: Albedo > 1.0: "
                  write(iulog,*) "SNICAR STATS: bnd_idx= ",bnd_idx
                  write (iulog,*) "SNICAR STATS: albout_lcl(bnd)= ",albout_lcl(bnd_idx), &
                     " albsfc_lcl(bnd_idx)= ",albsfc_lcl(bnd_idx)
                  write (iulog,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
                  write (iulog,*) "SNICAR STATS: coszen= ", coszen, " flg_slr= ", flg_slr_in

                  write (iulog,*) "SNICAR STATS: soot(-4)= ", mss_cnc_aer_lcl(-4,1)
                  write (iulog,*) "SNICAR STATS: soot(-3)= ", mss_cnc_aer_lcl(-3,1)
                  write (iulog,*) "SNICAR STATS: soot(-2)= ", mss_cnc_aer_lcl(-2,1)
                  write (iulog,*) "SNICAR STATS: soot(-1)= ", mss_cnc_aer_lcl(-1,1)
                  write (iulog,*) "SNICAR STATS: soot(0)= ", mss_cnc_aer_lcl(0,1)

                  write (iulog,*) "SNICAR STATS: L_snw(-4)= ", L_snw(-4)
                  write (iulog,*) "SNICAR STATS: L_snw(-3)= ", L_snw(-3)
                  write (iulog,*) "SNICAR STATS: L_snw(-2)= ", L_snw(-2)
                  write (iulog,*) "SNICAR STATS: L_snw(-1)= ", L_snw(-1)
                  write (iulog,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)

                  write (iulog,*) "SNICAR STATS: snw_rds(-4)= ", snw_rds(-4)
                  write (iulog,*) "SNICAR STATS: snw_rds(-3)= ", snw_rds(-3)
                  write (iulog,*) "SNICAR STATS: snw_rds(-2)= ", snw_rds(-2)
                  write (iulog,*) "SNICAR STATS: snw_rds(-1)= ", snw_rds(-1)
                  write (iulog,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(0)

                  CALL abort
               ENDIF
#endif
            ENDIF

         ENDDO   ! loop over wvl bands


         ! Weight output NIR albedo appropriately
         albout(1) = albout_lcl(1)
         flx_sum         = 0._r8
         DO bnd_idx= nir_bnd_bgn,nir_bnd_end
            flx_sum = flx_sum + flx_wgt(bnd_idx)*albout_lcl(bnd_idx)
         ENDDO
         albout(2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))

         ! Weight output NIR absorbed layer fluxes (flx_abs) appropriately
         flx_abs(:,1) = flx_abs_lcl(:,1)
         DO i=snl_top,1,1
            flx_sum = 0._r8
            DO bnd_idx= nir_bnd_bgn,nir_bnd_end
               flx_sum = flx_sum + flx_wgt(bnd_idx)*flx_abs_lcl(i,bnd_idx)
            ENDDO
            flx_abs(i,2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
         ENDDO

      ! If snow < minimum_snow, but > 0, and there is sun, set albedo to underlying surface albedo
      ELSEIF ( (coszen > 0._r8) .and. (h2osno_lcl < min_snw) .and. (h2osno_lcl > 0._r8) ) THEN
         albout(1) = albsfc(1)
         albout(2) = albsfc(2)

         ! There is either zero snow, or no sun
      ELSE
         albout(1) = 0._r8
         albout(2) = 0._r8
      ENDIF    ! IF column has snow and coszen > 0

   ! END associate

   END SUBROUTINE SNICAR_RT
   !-----------------------------------------------------------------------


   SUBROUTINE SNICAR_AD_RT (flg_snw_ice, flg_slr_in, &
                            coszen, snl, h2osno, frac_sno, &
                            h2osno_liq, h2osno_ice, snw_rds, &
                            mss_cnc_aer_in, albsfc, albout, flx_abs)
!-----------------------------------------------------------------------
!
! !DESCRIPTION:
!  Determine reflectance of, and vertically-resolved solar absorption in,
!  snow with impurities, with updated shortwave scheme
!
!  The multi-layer solution for multiple-scattering used here is from:
!  Briegleb, P. and Light, B.: A Delta-Eddington mutiple scattering
!  parameterization for solar radiation in the sea ice component of the
!  community climate system model, 2007.
!
!  The implementation of the SNICAR-AD model in ELM is described in:
!  Dang et al., Inter-comparison and improvement of 2-stream shortwave
!  radiative transfer models for unified treatment of cryospheric surfaces
!  in ESMs, in review, 2019
!
!  To USE this subtroutine, set use_snicar_ad = true in ELM
!
!  IF config_use_snicar_ad = true in MPAS-seaice
!  Snow on land and snow on sea ice will be treated
!  with the same model for their solar radiative properties.
!
!  The inputs and outputs are the same to SUBROUTINE SNICAR_RT
!
! !USES:
!-----------------------------------------------------------------------
! !ARGUMENTS:

   IMPLICIT NONE

   integer  , intent(in)  :: flg_snw_ice                  ! flag: =1 when called from CLM, =2 when called from CSIM
   integer  , intent(in)  :: flg_slr_in                   ! flag: =1 for direct-beam incident flux,=2 for diffuse incident flux
   real(r8) , intent(in)  :: coszen          ! cosine of solar zenith angle for next time step (col) [unitless]

   integer  , intent(in)  :: snl             !  negative number of snow layers (col) [nbr]
   real(r8) , intent(in)  :: h2osno          !  snow liquid water equivalent (col) [kg/m2]
   real(r8) , intent(in)  :: frac_sno        !  fraction of ground covered by snow (0 to 1)

   real(r8) , intent(in)  :: h2osno_liq     ( maxsnl+1:0 )   ! liquid water content (col,lyr) [kg/m2]
   real(r8) , intent(in)  :: h2osno_ice     ( maxsnl+1:0 )   ! ice content (col,lyr) [kg/m2]
   integer  , intent(in)  :: snw_rds        ( maxsnl+1:0 )   ! snow effective radius (col,lyr) [microns, m^-6]
   real(r8) , intent(in)  :: mss_cnc_aer_in ( maxsnl+1:0 , 1:sno_nbr_aer ) ! mass concentration of all aerosol species (col,lyr,aer) [kg/kg]
   real(r8) , intent(in)  :: albsfc         ( 1:numrad )     ! albedo of surface underlying snow (col,bnd) [frc]
   real(r8) , intent(out) :: albout         ( 1:numrad )     ! snow albedo, averaged into 2 bands (=0 IF no sun or no snow) (col,bnd) [frc]
   real(r8) , intent(out) :: flx_abs        ( maxsnl+1:1 , 1:numrad ) ! absorbed flux in each layer per unit flux incident (col, lyr, bnd)
   !
   ! !LOCAL VARIABLES:
   !
   ! variables for snow radiative transfer calculations

   ! Local variables representing single-column values of arrays:
   integer :: snl_lcl                                   ! negative number of snow layers [nbr]
   integer :: snw_rds_lcl(maxsnl+1:0)                   ! snow effective radius [m^-6]
   real(r8):: flx_slrd_lcl(1:numrad_snw)                ! direct beam incident irradiance [W/m2] (set to 1)
   real(r8):: flx_slri_lcl(1:numrad_snw)                ! diffuse incident irradiance [W/m2] (set to 1)
   real(r8):: mss_cnc_aer_lcl(maxsnl+1:0,1:sno_nbr_aer) ! aerosol mass concentration (lyr,aer_nbr) [kg/kg]
   real(r8):: h2osno_lcl                                ! total column snow mass [kg/m2]
   real(r8):: h2osno_liq_lcl(maxsnl+1:0)                ! liquid water mass [kg/m2]
   real(r8):: h2osno_ice_lcl(maxsnl+1:0)                ! ice mass [kg/m2]
   real(r8):: albsfc_lcl(1:numrad_snw)                  ! albedo of underlying surface [frc]
   real(r8):: ss_alb_snw_lcl(maxsnl+1:0)                ! single-scatter albedo of ice grains (lyr) [frc]
   real(r8):: asm_prm_snw_lcl(maxsnl+1:0)               ! asymmetry parameter of ice grains (lyr) [frc]
   real(r8):: ext_cff_mss_snw_lcl(maxsnl+1:0)           ! mass extinction coefficient of ice grains (lyr) [m2/kg]
   real(r8):: ss_alb_aer_lcl(sno_nbr_aer)               ! single-scatter albedo of aerosol species (aer_nbr) [frc]
   real(r8):: asm_prm_aer_lcl(sno_nbr_aer)              ! asymmetry parameter of aerosol species (aer_nbr) [frc]
   real(r8):: ext_cff_mss_aer_lcl(sno_nbr_aer)          ! mass extinction coefficient of aerosol species (aer_nbr) [m2/kg]

#ifdef MODAL_AER
   !mgf++
   real(r8) :: rds_bcint_lcl(maxsnl+1:0)       ! effective radius of within-ice BC [nm]
   real(r8) :: rds_bcext_lcl(maxsnl+1:0)       ! effective radius of external BC [nm]
   !mgf--
#endif


   ! Other local variables
   integer :: DELTA                               ! flag to USE Delta approximation (Joseph, 1976)
                                                  ! (1= USE, 0= don't USE)
   real(r8):: flx_wgt(1:numrad_snw)               ! weights applied to spectral bands,
                                                  ! specific to direct and diffuse cases (bnd) [frc]
   integer :: flg_nosnl                           ! flag: =1 IF there is snow, but zero snow layers,
                                                  ! =0 IF at least 1 snow layer [flg]
                                                  ! integer :: trip                                                                          ! flag: =1 to redo RT calculation IF result is unrealistic
                                                  ! integer :: flg_dover                                                                     ! defines conditions for RT redo (explained below)

   real(r8):: albedo                              ! temporary snow albedo [frc]
   real(r8):: flx_sum                             ! temporary summation variable for NIR weighting
   real(r8):: albout_lcl(numrad_snw)              ! snow albedo by band [frc]
   real(r8):: flx_abs_lcl(maxsnl+1:1,numrad_snw)  ! absorbed flux per unit incident flux at top of snowpack (lyr,bnd) [frc]

   real(r8):: L_snw(maxsnl+1:0)                   ! h2o mass (liquid+solid) in snow layer (lyr) [kg/m2]
   real(r8):: tau_snw(maxsnl+1:0)                 ! snow optical depth (lyr) [unitless]
   real(r8):: L_aer(maxsnl+1:0,sno_nbr_aer)       ! aerosol mass in snow layer (lyr,nbr_aer) [kg/m2]
   real(r8):: tau_aer(maxsnl+1:0,sno_nbr_aer)     ! aerosol optical depth (lyr,nbr_aer) [unitless]
   real(r8):: tau_sum                             ! cumulative (snow+aerosol) optical depth [unitless]
   real(r8):: tau_elm(maxsnl+1:0)                 ! column optical depth from layer bottom to snowpack top (lyr) [unitless]
   real(r8):: omega_sum                           ! temporary summation of single-scatter albedo of all aerosols [frc]
   real(r8):: g_sum                               ! temporary summation of asymmetry parameter of all aerosols [frc]

   real(r8):: tau(maxsnl+1:0)                     ! weighted optical depth of snow+aerosol layer (lyr) [unitless]
   real(r8):: omega(maxsnl+1:0)                   ! weighted single-scatter albedo of snow+aerosol layer (lyr) [frc]
   real(r8):: g(maxsnl+1:0)                       ! weighted asymmetry parameter of snow+aerosol layer (lyr) [frc]
   real(r8):: tau_star(maxsnl+1:0)                ! transformed (i.e. Delta-Eddington) optical depth of snow+aerosol layer
                                                  ! (lyr) [unitless]
   real(r8):: omega_star(maxsnl+1:0)              ! transformed (i.e. Delta-Eddington) SSA of snow+aerosol layer (lyr) [frc]
   real(r8):: g_star(maxsnl+1:0)                  ! transformed (i.e. Delta-Eddington) asymmetry paramater of snow+aerosol layer
                                                  ! (lyr) [frc]

                                                  ! integer :: c_idx                                                                         ! column indices [idx]
   integer :: bnd_idx                             ! spectral band index (1 <= bnd_idx <= numrad_snw) [idx]
   integer :: rds_idx                             ! snow effective radius index for retrieving
                                                  ! Mie parameters from lookup table [idx]
   integer :: snl_btm                             ! index of bottom snow layer (0) [idx]
   integer :: snl_top                             ! index of top snow layer (-4 to 0) [idx]
   integer :: fc                                  ! column filter index
   integer :: i                                   ! layer index [idx]
   integer :: j                                   ! aerosol number index [idx]
   integer :: m                                   ! secondary layer index [idx]
   integer :: nint_snw_rds_min                    ! nearest integer value of snw_rds_min

   real(r8):: F_abs(maxsnl+1:0)                   ! net absorbed radiative energy (lyr) [W/m^2]
   real(r8):: F_abs_sum                           ! total absorbed energy in column [W/m^2]
   real(r8):: F_sfc_pls                           ! upward radiative flux at snowpack top [W/m^2]
   real(r8):: F_btm_net                           ! net flux at bottom of snowpack [W/m^2]
   real(r8):: energy_sum                          ! sum of all energy terms; should be 0.0 [W/m^2]
   real(r8):: mu_not                              ! cosine of solar zenith angle (used locally) [frc]

   integer :: err_idx                             ! counter for number of times through error loop [nbr]
   real(r8):: pi                                  ! 3.1415...

   integer :: snw_shp_lcl(maxsnl+1:0)             ! Snow grain shape option:
                                                  ! 1=sphere; 2=spheroid; 3=hexagonal plate; 4=koch snowflake
   real(r8):: snw_fs_lcl(maxsnl+1:0)              ! Shape factor: ratio of nonspherical grain effective radii to that of equal-volume sphere
                                                  ! 0=USE recommended default value
                                                  ! others(0<fs<1)= USE user-specified value
                                                  ! only activated when sno_shp > 1 (i.e. nonspherical)
   real(r8):: snw_ar_lcl(maxsnl+1:0)              ! % Aspect ratio: ratio of grain width to length
                                                  ! 0=USE recommended default value
                                                  ! others(0.1<fs<20)= USE user-specified value
                                                  ! only activated when sno_shp > 1 (i.e. nonspherical)
   real(r8):: &
       diam_ice           , & ! effective snow grain diameter
       fs_sphd            , & ! shape factor for spheroid
       fs_hex0            , & ! shape factor for hexagonal plate
       fs_hex             , & ! shape factor for hexagonal plate (reference)
       fs_koch            , & ! shape factor for koch snowflake
       AR_tmp             , & ! aspect ratio for spheroid
       g_ice_Cg_tmp(7)    , & ! temporary for calculation of asymetry factor
       gg_ice_F07_tmp(7)  , & ! temporary for calculation of asymetry factor
       g_ice_F07          , & ! temporary for calculation of asymetry factor
       g_ice              , & ! asymmetry factor
       gg_F07_intp        , & ! temporary for calculation of asymetry factor (interpolated)
       g_Cg_intp          , & ! temporary for calculation of asymetry factor  (interpolated)
       R_1_omega_tmp      , & ! temporary for dust-snow mixing calculation
       C_dust_total           ! dust concentration

   integer :: atm_type_index  ! index for atmospheric type
   integer :: slr_zen         ! integer value of solar zenith angle

   ! SNICAR_AD new variables, follow sea-ice shortwave conventions
   real(r8):: &
      trndir(maxsnl+1:1)  , & ! solar beam down transmission from top
      trntdr(maxsnl+1:1)  , & ! total transmission to direct beam for layers above
      trndif(maxsnl+1:1)  , & ! diffuse transmission to diffuse beam for layers above
      rupdir(maxsnl+1:1)  , & ! reflectivity to direct radiation for layers below
      rupdif(maxsnl+1:1)  , & ! reflectivity to diffuse radiation for layers below
      rdndif(maxsnl+1:1)  , & ! reflectivity to diffuse radiation for layers above
      dfdir(maxsnl+1:1)   , & ! down-up flux at interface due to direct beam at top surface
      dfdif(maxsnl+1:1)   , & ! down-up flux at interface due to diffuse beam at top surface
      dftmp(maxsnl+1:1)       ! temporary variable for down-up flux at interface

   real(r8):: &
      rdir(maxsnl+1:0)    , & ! layer reflectivity to direct radiation
      rdif_a(maxsnl+1:0)  , & ! layer reflectivity to diffuse radiation from above
      rdif_b(maxsnl+1:0)  , & ! layer reflectivity to diffuse radiation from below
      tdir(maxsnl+1:0)    , & ! layer transmission to direct radiation (solar beam + diffuse)
      tdif_a(maxsnl+1:0)  , & ! layer transmission to diffuse radiation from above
      tdif_b(maxsnl+1:0)  , & ! layer transmission to diffuse radiation from below
      trnlay(maxsnl+1:0)      ! solar beam transm for layer (direct beam only)

   real(r8):: &
       ts       , & ! layer delta-scaled extinction optical depth
       ws       , & ! layer delta-scaled single scattering albedo
       gs       , & ! layer delta-scaled asymmetry parameter
       extins   , & ! extinction
       alp      , & ! temporary for alpha
       gam      , & ! temporary for agamm
       amg      , & ! alp - gam
       apg      , & ! alp + gam
       ue       , & ! temporary for u
       refk     , & ! interface multiple scattering
       refkp1   , & ! interface multiple scattering for k+1
       refkm1   , & ! interface multiple scattering for k-1
       tdrrdir  , & ! direct tran times layer direct ref
       tdndif       ! total down diffuse = tot tran - direct tran

   real(r8) :: &
       alpha    , & ! term in direct reflectivity and transmissivity
       agamm    , & ! term in direct reflectivity and transmissivity
       el       , & ! term in alpha,agamm,n,u
       taus     , & ! scaled extinction optical depth
       omgs     , & ! scaled single particle scattering albedo
       asys     , & ! scaled asymmetry parameter
       u        , & ! term in diffuse reflectivity and transmissivity
       n        , & ! term in diffuse reflectivity and transmissivity
       lm       , & ! temporary for el
       mu       , & ! cosine solar zenith for either snow or water
       ne           ! temporary for n

   ! perpendicular and parallel relative to plane of incidence and scattering
   real(r8) :: &
       R1       , & ! perpendicular polarization reflection amplitude
       R2       , & ! parallel polarization reflection amplitude
       T1       , & ! perpendicular polarization transmission amplitude
       T2       , & ! parallel polarization transmission amplitude
       Rf_dir_a , & ! fresnel reflection to direct radiation
       Tf_dir_a , & ! fresnel transmission to direct radiation
       Rf_dif_a , & ! fresnel reflection to diff radiation from above
       Rf_dif_b , & ! fresnel reflection to diff radiation from below
       Tf_dif_a , & ! fresnel transmission to diff radiation from above
       Tf_dif_b     ! fresnel transmission to diff radiation from below

   real(r8) :: &
       gwt      , & ! gaussian weight
       swt      , & ! sum of weights
       trn      , & ! layer transmission
       rdr      , & ! rdir for gaussian integration
       tdr      , & ! tdir for gaussian integration
       smr      , & ! accumulator for rdif gaussian integration
       smt      , & ! accumulator for tdif gaussian integration
       exp_min      ! minimum exponential value

   integer :: &
       ng             , & ! gaussian integration index
       snl_btm_itf    , & ! index of bottom snow layer interfaces (1) [idx]
       ngmax = 8          ! gaussian integration index

   ! Gaussian integration angle and coefficients
   real(r8) :: difgauspt(1:8) , difgauswt(1:8)

   ! constants used in algorithm
   real(r8) :: &
       c0      = 0.0_r8     , &
       c1      = 1.0_r8     , &
       c3      = 3.0_r8     , &
       c4      = 4.0_r8     , &
       c6      = 6.0_r8     , &
       cp01    = 0.01_r8    , &
       cp5     = 0.5_r8     , &
       cp75    = 0.75_r8    , &
       c1p5    = 1.5_r8     , &
       trmin   = 0.001_r8   , &
       argmax  = 10.0_r8       ! maximum argument of exponential

   ! cconstant coefficients used for SZA parameterization
   real(r8) :: &
       sza_a0 =  0.085730_r8 , &
       sza_a1 = -0.630883_r8 , &
       sza_a2 =  1.303723_r8 , &
       sza_b0 =  1.467291_r8 , &
       sza_b1 = -3.338043_r8 , &
       sza_b2 =  6.807489_r8 , &
       puny   =  1.0e-11_r8  , &
       mu_75  =  0.2588_r8     ! cosine of 75 degree

   ! coefficients used for SZA parameterization
   real(r8) :: &
       sza_c1          , & ! coefficient, SZA parameteirzation
       sza_c0          , & ! coefficient, SZA parameterization
       sza_factor      , & ! factor used to adjust NIR direct albedo
       flx_sza_adjust  , & ! direct NIR flux adjustment from sza_factor
       mu0                 ! incident solar zenith angle

     !-----------------------------------------------------------------------
#ifdef MODAL_AER
   !mgf++
   integer :: idx_bcint_icerds  ! index of ice effective radius for optical properties lookup table
   integer :: idx_bcint_nclrds  ! index of within-ice BC effective radius for optical properties lookup table
   integer :: idx_bcext_nclrds  ! index of external BC effective radius for optical properties lookup table
   real(r8):: enh_fct           ! extinction/absorption enhancement factor for within-ice BC
   real(r8):: tmp1              ! temporary variable
   !mgf--
#endif

   ! Constants for non-spherical ice particles and dust-snow internal mixing
   real(r8) :: g_b2(7)
   real(r8) :: g_b1(7)
   real(r8) :: g_b0(7)
   real(r8) :: g_F07_c2(7)
   real(r8) :: g_F07_c1(7)
   real(r8) :: g_F07_c0(7)
   real(r8) :: g_F07_p2(7)
   real(r8) :: g_F07_p1(7)
   real(r8) :: g_F07_p0(7)
   real(r8) :: dust_clear_d0(3)
   real(r8) :: dust_clear_d1(3)
   real(r8) :: dust_clear_d2(3)
   real(r8) :: dust_cloudy_d0(3)
   real(r8) :: dust_cloudy_d1(3)
   real(r8) :: dust_cloudy_d2(3)

   !!! factors for considering snow grain shape
   data g_b0(:) / 9.76029E-01_r8, 9.67798E-01_r8, 1.00111E+00_r8, 1.00224E+00_r8,&
                  9.64295E-01_r8, 9.97475E-01_r8, 9.97475E-01_r8/
   data g_b1(:) / 5.21042E-01_r8, 4.96181E-01_r8, 1.83711E-01_r8, 1.37082E-01_r8,&
                  5.50598E-02_r8, 8.48743E-02_r8, 8.48743E-02_r8/
   data g_b2(:) /-2.66792E-04_r8, 1.14088E-03_r8, 2.37011E-04_r8,-2.35905E-04_r8,&
                  8.40449E-04_r8,-4.71484E-04_r8,-4.71484E-04_r8/

   data g_F07_c2(:) / 1.349959E-1_r8, 1.115697E-1_r8, 9.853958E-2_r8, 5.557793E-2_r8,&
                     -1.233493E-1_r8, 0.0_r8, 0.0_r8/
   data g_F07_c1(:) /-3.987320E-1_r8,-3.723287E-1_r8,-3.924784E-1_r8,-3.259404E-1_r8,&
                      4.429054E-2_r8,-1.726586E-1_r8,-1.726586E-1_r8/
   data g_F07_c0(:) / 7.938904E-1_r8, 8.030084E-1_r8, 8.513932E-1_r8, 8.692241E-1_r8,&
                      7.085850E-1_r8, 6.412701E-1_r8, 6.412701E-1_r8/
   data g_F07_p2(:) / 3.165543E-3_r8, 2.014810E-3_r8, 1.780838E-3_r8, 6.987734E-4_r8,&
                     -1.882932E-2_r8,-2.277872E-2_r8,-2.277872E-2_r8/
   data g_F07_p1(:) / 1.140557E-1_r8, 1.143152E-1_r8, 1.143814E-1_r8, 1.071238E-1_r8,&
                      1.353873E-1_r8, 1.914431E-1_r8, 1.914431E-1_r8/
   data g_F07_p0(:) / 5.292852E-1_r8, 5.425909E-1_r8, 5.601598E-1_r8, 6.023407E-1_r8,&
                      6.473899E-1_r8, 4.634944E-1_r8, 4.634944E-1_r8/

   !!! factors for considring dust-snow internal mixing
   data dust_clear_d0(:) /1.0413E+00_r8,1.0168E+00_r8,1.0189E+00_r8/
   data dust_clear_d1(:) /1.0016E+00_r8,1.0070E+00_r8,1.0840E+00_r8/
   data dust_clear_d2(:) /2.4208E-01_r8,1.5300E-03_r8,1.1230E-04_r8/

   data dust_cloudy_d0(:) /1.0388E+00_r8,1.0167E+00_r8,1.0189E+00_r8/
   data dust_cloudy_d1(:) /1.0015E+00_r8,1.0061E+00_r8,1.0823E+00_r8/
   data dust_cloudy_d2(:) /2.5973E-01_r8,1.6200E-03_r8,1.1721E-04_r8/

   ! Enforce expected array sizes

  !  associate(&
  !       snl      => col_pp%snl           , & ! Input:  [integer (:)]  negative number of snow layers (col) [nbr]
  !       h2osno   => col_ws%h2osno        , & ! Input:  [real(r8) (:)]  snow liquid water equivalent (col) [kg/m2]
  !       frac_sno => col_ws%frac_sno_eff    & ! Input:  [real(r8) (:)]  fraction of ground covered by snow (0 to 1)
  !       )

      ! Define constants
      pi = SHR_CONST_PI
      nint_snw_rds_min = nint(snw_rds_min)

      ! always USE Delta approximation for snow
      DELTA = 1

      !Gaussian integration angle and coefficients for diffuse radiation
      difgauspt(1:8)     & ! gaussian angles (radians)
        = (/ 0.9894009_r8,  0.9445750_r8, &
             0.8656312_r8,  0.7554044_r8, &
             0.6178762_r8,  0.4580168_r8, &
             0.2816036_r8,  0.0950125_r8/)
      difgauswt(1:8)     & ! gaussian weights
        = (/ 0.0271525_r8,  0.0622535_r8, &
             0.0951585_r8,  0.1246290_r8, &
             0.1495960_r8,  0.1691565_r8, &
             0.1826034_r8,  0.1894506_r8/)

      snw_shp_lcl(:) = snow_shape_sphere
      snw_fs_lcl(:)  = 0._r8
      snw_ar_lcl(:)  = 0._r8
      atm_type_index = atm_type_default

      ! Define snow grain shape
      IF (trim(snow_shape) == 'sphere') THEN
        snw_shp_lcl(:) = snow_shape_sphere
      ELSEIF (trim(snow_shape) == 'spheroid') THEN
        snw_shp_lcl(:) = snow_shape_spheroid
      ELSEIF (trim(snow_shape) == 'hexagonal_plate') THEN
        snw_shp_lcl(:) = snow_shape_hexagonal_plate
      ELSEIF (trim(snow_shape) == 'koch_snowflake') THEN
        snw_shp_lcl(:) = snow_shape_koch_snowflake
      ELSE
         IF (p_is_master) THEN
            write(iulog,*) "snow_shape = ", snow_shape
            CALL abort
         ENDIF
      ENDIF

      ! Define atmospheric type
      IF (trim(snicar_atm_type) == 'default') THEN
        atm_type_index = atm_type_default
      ELSEIF (trim(snicar_atm_type) == 'mid-latitude_winter') THEN
        atm_type_index = atm_type_mid_latitude_winter
      ELSEIF (trim(snicar_atm_type) == 'mid-latitude_summer') THEN
        atm_type_index = atm_type_mid_latitude_summer
      ELSEIF (trim(snicar_atm_type) == 'sub-Arctic_winter') THEN
        atm_type_index = atm_type_sub_Arctic_winter
      ELSEIF (trim(snicar_atm_type) == 'sub-Arctic_summer') THEN
        atm_type_index = atm_type_sub_Arctic_summer
      ELSEIF (trim(snicar_atm_type) == 'summit_Greenland') THEN
        atm_type_index = atm_type_summit_Greenland
      ELSEIF (trim(snicar_atm_type) == 'high_mountain') THEN
        atm_type_index = atm_type_high_mountain
      ELSE
         IF (p_is_master) THEN
            write(iulog,*) "snicar_atm_type = ", snicar_atm_type
            CALL abort
         ENDIF
      ENDIF

      ! (when called from CSIM, there is only one column)

      ! Zero absorbed radiative fluxes:
      DO i=maxsnl+1,1,1
         flx_abs_lcl(:,:)   = 0._r8
         flx_abs(i,:) = 0._r8
      ENDDO

      ! set snow/ice mass to be used for RT:
      IF (flg_snw_ice == 1) THEN
         h2osno_lcl = h2osno
      ELSE
         h2osno_lcl = h2osno_ice(0)
      ENDIF

      ! Qualifier for computing snow RT:
      !  1) sunlight from atmosphere model
      !  2) minimum amount of snow on ground.
      !     Otherwise, set snow albedo to zero
      IF ((coszen > 0._r8) .and. (h2osno_lcl > min_snw) ) THEN

         ! Set variables specific to ELM
         IF (flg_snw_ice == 1) THEN
            ! If there is snow, but zero snow layers, we must create a layer locally.
            ! This layer is presumed to have the fresh snow effective radius.
            IF (snl > -1) THEN
               flg_nosnl         =  1
               snl_lcl           =  -1
               h2osno_ice_lcl(0) =  h2osno_lcl
               h2osno_liq_lcl(0) =  0._r8
               snw_rds_lcl(0)    =  nint_snw_rds_min
            ELSE
               flg_nosnl         =  0
               snl_lcl           =  snl
               h2osno_liq_lcl(:) =  h2osno_liq(:)
               h2osno_ice_lcl(:) =  h2osno_ice(:)
               snw_rds_lcl(:)    =  snw_rds(:)
            ENDIF

            snl_btm   = 0
            snl_top   = snl_lcl+1

            ! Set variables specific to CSIM
         ELSE
            flg_nosnl         = 0
            snl_lcl           = -1
            h2osno_liq_lcl(:) = h2osno_liq(:)
            h2osno_ice_lcl(:) = h2osno_ice(:)
            snw_rds_lcl(:)    = snw_rds(:)
            snl_btm           = 0
            snl_top           = 0
         ENDIF ! END IF flg_snw_ice == 1

#ifdef MODAL_AER
         !mgf++
         !
         ! Assume fixed BC effective radii of 100nm. This is close to
         ! the effective radius of 95nm (number median radius of
         ! 40nm) assumed for freshly-emitted BC in MAM.  Future
         ! implementations may prognose the BC effective radius in
         ! snow.
         rds_bcint_lcl(:)  =  100._r8
         rds_bcext_lcl(:)  =  100._r8
         !mgf--
#endif

         ! Set local aerosol array
         DO j=1,sno_nbr_aer
            mss_cnc_aer_lcl(:,j) = mss_cnc_aer_in(:,j)
         ENDDO

         ! Set spectral underlying surface albedos to their corresponding VIS or NIR albedos
         albsfc_lcl(1)                       = albsfc(1)
         albsfc_lcl(nir_bnd_bgn:nir_bnd_end) = albsfc(2)

         ! Error check for snow grain size:
         IF (p_is_master) THEN
            DO i=snl_top,snl_btm,1
               IF ((snw_rds_lcl(i) < snw_rds_min_tbl) .or. (snw_rds_lcl(i) > snw_rds_max_tbl)) THEN
                  write (iulog,*) "SNICAR ERROR: snow grain radius of ", snw_rds_lcl(i), " out of bounds."
                  write (iulog,*) "flg_snw_ice= ", flg_snw_ice
                  write (iulog,*) " level: ", i, " snl(c)= ", snl_lcl
                  write (iulog,*) "h2osno(c)= ", h2osno_lcl
                  CALL abort
               ENDIF
            ENDDO
         ENDIF

         ! Incident flux weighting parameters
         !  - sum of all VIS bands must equal 1
         !  - sum of all NIR bands must equal 1
         !
         ! Spectral bands (5-band CASE)
         !  Band 1: 0.3-0.7um (VIS)
         !  Band 2: 0.7-1.0um (NIR)
         !  Band 3: 1.0-1.2um (NIR)
         !  Band 4: 1.2-1.5um (NIR)
         !  Band 5: 1.5-5.0um (NIR)
         !
         ! The following weights are appropriate for surface-incident flux in a mid-latitude winter atmosphere
         !
         ! 3-band weights
         IF (numrad_snw==3) THEN
            ! Direct:
            IF (flg_slr_in == 1) THEN
               flx_wgt(1) = 1._r8
               flx_wgt(2) = 0.66628670195247_r8
               flx_wgt(3) = 0.33371329804753_r8
               ! Diffuse:
            ELSEIF (flg_slr_in == 2) THEN
               flx_wgt(1) = 1._r8
               flx_wgt(2) = 0.77887652162877_r8
               flx_wgt(3) = 0.22112347837123_r8
            ENDIF

            ! 5-band weights
         ELSEIF(numrad_snw==5) THEN
            ! Direct:
            IF (flg_slr_in == 1) THEN
              IF (atm_type_index == atm_type_default) THEN
                 flx_wgt(1) = 1._r8
                 flx_wgt(2) = 0.49352158521175_r8
                 flx_wgt(3) = 0.18099494230665_r8
                 flx_wgt(4) = 0.12094898498813_r8
                 flx_wgt(5) = 0.20453448749347_r8
              ELSE
                 slr_zen = nint(acos(coszen) * 180._r8 / pi)
                 IF (slr_zen>89) THEN
                    slr_zen = 89
                 ENDIF
                 flx_wgt(1) = 1._r8
                 flx_wgt(2) = flx_wgt_dir(atm_type_index, slr_zen+1, 2)
                 flx_wgt(3) = flx_wgt_dir(atm_type_index, slr_zen+1, 3)
                 flx_wgt(4) = flx_wgt_dir(atm_type_index, slr_zen+1, 4)
                 flx_wgt(5) = flx_wgt_dir(atm_type_index, slr_zen+1, 5)
              ENDIF

               ! Diffuse:
            ELSEIF (flg_slr_in == 2) THEN
               IF  (atm_type_index == atm_type_default) THEN
                 flx_wgt(1) = 1._r8
                 flx_wgt(2) = 0.58581507618433_r8
                 flx_wgt(3) = 0.20156903770812_r8
                 flx_wgt(4) = 0.10917889346386_r8
                 flx_wgt(5) = 0.10343699264369_r8
              ELSE
                 flx_wgt(1) = 1._r8
                 flx_wgt(2) = flx_wgt_dif(atm_type_index, 2)
                 flx_wgt(3) = flx_wgt_dif(atm_type_index, 3)
                 flx_wgt(4) = flx_wgt_dif(atm_type_index, 4)
                 flx_wgt(5) = flx_wgt_dif(atm_type_index, 5)
              ENDIF
            ENDIF
         ENDIF ! END IF numrad_snw

         ! Loop over snow spectral bands

         exp_min = exp(-argmax)
         DO bnd_idx = 1,numrad_snw

           ! note that we can remove flg_dover since this algorithm is
           ! stable for mu_not > 0.01

           ! mu_not is cosine solar zenith angle above the fresnel level; make
           ! sure mu_not is large enough for stable and meaningful radiation
           ! solution: .01 is like sun just touching horizon with its lower edge
           ! equivalent to mu0 in sea-ice shortwave model ice_shortwave.F90
            mu_not = max(coszen, cp01)


            ! Set direct or diffuse incident irradiance to 1
            ! (This has to be within the bnd loop because mu_not is adjusted in rare cases)
            IF (flg_slr_in == 1) THEN
               flx_slrd_lcl(bnd_idx) = 1._r8/(mu_not*pi) ! this corresponds to incident irradiance of 1.0
               flx_slri_lcl(bnd_idx) = 0._r8
            ELSE
               flx_slrd_lcl(bnd_idx) = 0._r8
               flx_slri_lcl(bnd_idx) = 1._r8
            ENDIF

            ! Pre-emptive error handling: aerosols can reap havoc on these absorptive bands.
            ! Since extremely high soot concentrations have a negligible effect on these bands, zero them.
            IF ( (numrad_snw == 5).and.((bnd_idx == 5).or.(bnd_idx == 4)) ) THEN
               mss_cnc_aer_lcl(:,:) = 0._r8
            ENDIF

            IF ( (numrad_snw == 3).and.(bnd_idx == 3) ) THEN
               mss_cnc_aer_lcl(:,:) = 0._r8
            ENDIF

            ! Define local Mie parameters based on snow grain size and aerosol species,
            !  retrieved from a lookup table.
            IF (flg_slr_in == 1) THEN
               DO i=snl_top,snl_btm,1
                  rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
                  ! snow optical properties (direct radiation)
                  ss_alb_snw_lcl(i)      = ss_alb_snw_drc(rds_idx,bnd_idx)
                  asm_prm_snw_lcl(i)     = asm_prm_snw_drc(rds_idx,bnd_idx)
                  ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_drc(rds_idx,bnd_idx)
               ENDDO
            ELSEIF (flg_slr_in == 2) THEN
               DO i=snl_top,snl_btm,1
                  rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
                  ! snow optical properties (diffuse radiation)
                  ss_alb_snw_lcl(i)      = ss_alb_snw_dfs(rds_idx,bnd_idx)
                  asm_prm_snw_lcl(i)     = asm_prm_snw_dfs(rds_idx,bnd_idx)
                  ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_dfs(rds_idx,bnd_idx)
               ENDDO
            ENDIF

           ! Calculate the asymetry factors under different snow grain shapes
            DO i=snl_top,snl_btm,1
               IF(snw_shp_lcl(i) == snow_shape_spheroid) THEN ! spheroid
                  diam_ice = 2._r8*snw_rds_lcl(i)
                  IF(snw_fs_lcl(i) == 0._r8) THEN
                     fs_sphd = 0.929_r8
                  ELSE
                     fs_sphd = snw_fs_lcl(i)
                  ENDIF
                  fs_hex = 0.788_r8
                  IF(snw_ar_lcl(i) == 0._r8) THEN
                     AR_tmp = 0.5_r8
                  ELSE
                     AR_tmp = snw_ar_lcl(i)
                  ENDIF
                  g_ice_Cg_tmp = g_b0 * ((fs_sphd/fs_hex)**g_b1) * (diam_ice**g_b2)
                  gg_ice_F07_tmp = g_F07_c0 + g_F07_c1 * AR_tmp + g_F07_c2 * (AR_tmp**2)
               ELSEIF(snw_shp_lcl(i) == snow_shape_hexagonal_plate) THEN ! hexagonal plate
                  diam_ice = 2._r8*snw_rds_lcl(i)
                  IF(snw_fs_lcl(i) == 0._r8) THEN
                     fs_hex0 = 0.788_r8
                  ELSE
                     fs_hex0 = snw_fs_lcl(i)
                  ENDIF
                  fs_hex = 0.788_r8
                  IF(snw_ar_lcl(i) == 0._r8) THEN
                    AR_tmp = 2.5_r8
                  ELSE
                    AR_tmp = snw_ar_lcl(i)
                  ENDIF
                  g_ice_Cg_tmp = g_b0 * ((fs_hex0/fs_hex)**g_b1) * (diam_ice**g_b2)
                  gg_ice_F07_tmp = g_F07_p0 + g_F07_p1 * log(AR_tmp) + g_F07_p2 * ((log(AR_tmp))**2)
               ELSEIF(snw_shp_lcl(i) == snow_shape_koch_snowflake) THEN ! Koch snowflake
                  diam_ice = 2._r8 * snw_rds_lcl(i) /0.544_r8
                  IF(snw_fs_lcl(i) == 0._r8) THEN
                     fs_koch = 0.712_r8
                  ELSE
                     fs_koch = snw_fs_lcl(i)
                  ENDIF
                  fs_hex = 0.788_r8
                  IF(snw_ar_lcl(i) == 0._r8) THEN
                     AR_tmp = 2.5_r8
                  ELSE
                     AR_tmp = snw_ar_lcl(i)
                  ENDIF
                  g_ice_Cg_tmp = g_b0 * ((fs_koch/fs_hex)**g_b1) * (diam_ice**g_b2)
                  gg_ice_F07_tmp = g_F07_p0 + g_F07_p1 * log(AR_tmp) + g_F07_p2 * ((log(AR_tmp))**2)
              ENDIF

              ! Linear interpolation for calculating the asymetry factor at band_idx.
              IF(snw_shp_lcl(i) > 1) THEN
                IF(bnd_idx == 1) THEN
                  g_Cg_intp   = (g_ice_Cg_tmp(2)-g_ice_Cg_tmp(1))/(1.055_r8-0.475_r8)*(0.5_r8-0.475_r8) +g_ice_Cg_tmp(1)
                  gg_F07_intp = (gg_ice_F07_tmp(2)-gg_ice_F07_tmp(1))/(1.055_r8-0.475_r8)*(0.5_r8-0.475_r8)+gg_ice_F07_tmp(1)
                ELSEIF(bnd_idx == 2) THEN
                  g_Cg_intp   = (g_ice_Cg_tmp(2)-g_ice_Cg_tmp(1))/(1.055_r8-0.475_r8)*(0.85_r8-0.475_r8)+g_ice_Cg_tmp(1)
                gg_F07_intp = (gg_ice_F07_tmp(2)-gg_ice_F07_tmp(1))/(1.055_r8-0.475_r8)*(0.85_r8-0.475_r8)+gg_ice_F07_tmp(1)
                ELSEIF(bnd_idx == 3) THEN
                  g_Cg_intp   = (g_ice_Cg_tmp(3)-g_ice_Cg_tmp(2))/(1.655_r8-1.055_r8)*(1.1_r8-1.055_r8)&
                                +g_ice_Cg_tmp(2)
                  gg_F07_intp = (gg_ice_F07_tmp(3)-gg_ice_F07_tmp(2))/(1.655_r8-1.055_r8)*(1.1_r8-1.055_r8)&
                                +gg_ice_F07_tmp(2)
                ELSEIF(bnd_idx == 4) THEN
                  g_Cg_intp   = (g_ice_Cg_tmp(3)-g_ice_Cg_tmp(2))/(1.655_r8-1.055_r8)*(1.35_r8-1.055_r8)&
                                +g_ice_Cg_tmp(2)
                  gg_F07_intp = (gg_ice_F07_tmp(3)-gg_ice_F07_tmp(2))/(1.655_r8-1.055_r8)*(1.35_r8-1.055_r8)&
                                +gg_ice_F07_tmp(2)
                ELSEIF(bnd_idx == 5) THEN
                  g_Cg_intp   = (g_ice_Cg_tmp(6)-g_ice_Cg_tmp(5))/(3.75_r8-3.0_r8)*(3.25_r8-3.0_r8)&
                                +g_ice_Cg_tmp(5)
                  gg_F07_intp = (gg_ice_F07_tmp(6)-gg_ice_F07_tmp(5))/(3.75_r8-3.0_r8)*(3.25_r8-3.0_r8)&
                                +gg_ice_F07_tmp(5)
                ENDIF
                g_ice_F07 = gg_F07_intp + (1._r8 - gg_F07_intp) / ss_alb_snw_lcl(i) / 2._r8
                g_ice = g_ice_F07 * g_Cg_intp
                asm_prm_snw_lcl(i) = g_ice
              ENDIF

              IF(asm_prm_snw_lcl(i) > 0.99_r8) THEN
                asm_prm_snw_lcl(i) = 0.99_r8
              ENDIF

            ENDDO
            !!!-END

   !H. Wang
            ! aerosol species 1 optical properties
            ! ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx)
            ! asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx)
            ! ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx)

            ! aerosol species 2 optical properties
            ! ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx)
            ! asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx)
            ! ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx)
   !H. Wang
            ! aerosol species 3 optical properties
            ss_alb_aer_lcl(3)        = ss_alb_oc1(bnd_idx)
            asm_prm_aer_lcl(3)       = asm_prm_oc1(bnd_idx)
            ext_cff_mss_aer_lcl(3)   = ext_cff_mss_oc1(bnd_idx)

            ! aerosol species 4 optical properties
            ss_alb_aer_lcl(4)        = ss_alb_oc2(bnd_idx)
            asm_prm_aer_lcl(4)       = asm_prm_oc2(bnd_idx)
            ext_cff_mss_aer_lcl(4)   = ext_cff_mss_oc2(bnd_idx)

            ! aerosol species 5 optical properties
            ss_alb_aer_lcl(5)        = ss_alb_dst1(bnd_idx)
            asm_prm_aer_lcl(5)       = asm_prm_dst1(bnd_idx)
            ext_cff_mss_aer_lcl(5)   = ext_cff_mss_dst1(bnd_idx)

            ! aerosol species 6 optical properties
            ss_alb_aer_lcl(6)        = ss_alb_dst2(bnd_idx)
            asm_prm_aer_lcl(6)       = asm_prm_dst2(bnd_idx)
            ext_cff_mss_aer_lcl(6)   = ext_cff_mss_dst2(bnd_idx)

            ! aerosol species 7 optical properties
            ss_alb_aer_lcl(7)        = ss_alb_dst3(bnd_idx)
            asm_prm_aer_lcl(7)       = asm_prm_dst3(bnd_idx)
            ext_cff_mss_aer_lcl(7)   = ext_cff_mss_dst3(bnd_idx)

            ! aerosol species 8 optical properties
            ss_alb_aer_lcl(8)        = ss_alb_dst4(bnd_idx)
            asm_prm_aer_lcl(8)       = asm_prm_dst4(bnd_idx)
            ext_cff_mss_aer_lcl(8)   = ext_cff_mss_dst4(bnd_idx)


            ! 1. snow and aerosol layer column mass (L_snw, L_aer [kg/m^2])
            ! 2. optical Depths (tau_snw, tau_aer)
            ! 3. weighted Mie properties (tau, omega, g)

            ! Weighted Mie parameters of each layer
            DO i=snl_top,snl_btm,1
#ifdef MODAL_AER
               !mgf++ within-ice and external BC optical properties
               !
               ! Lookup table indices for BC optical properties,
               ! dependent on snow grain size and BC particle
               ! size.

               ! valid for 25 < snw_rds < 1625 um:
               IF (snw_rds_lcl(i) < 125) THEN
                  tmp1 = snw_rds_lcl(i)/50
                  idx_bcint_icerds = nint(tmp1)
               ELSEIF (snw_rds_lcl(i) < 175) THEN
                  idx_bcint_icerds = 2
               ELSE
                  tmp1 = (snw_rds_lcl(i)/250)+2
                  idx_bcint_icerds = nint(tmp1)
               ENDIF

               ! valid for 25 < bc_rds < 525 nm
               idx_bcint_nclrds = nint(rds_bcint_lcl(i)/50)
               idx_bcext_nclrds = nint(rds_bcext_lcl(i)/50)

               ! check bounds:
               IF (idx_bcint_icerds < idx_bcint_icerds_min) idx_bcint_icerds = idx_bcint_icerds_min
               IF (idx_bcint_icerds > idx_bcint_icerds_max) idx_bcint_icerds = idx_bcint_icerds_max
               IF (idx_bcint_nclrds < idx_bc_nclrds_min) idx_bcint_nclrds = idx_bc_nclrds_min
               IF (idx_bcint_nclrds > idx_bc_nclrds_max) idx_bcint_nclrds = idx_bc_nclrds_max
               IF (idx_bcext_nclrds < idx_bc_nclrds_min) idx_bcext_nclrds = idx_bc_nclrds_min
               IF (idx_bcext_nclrds > idx_bc_nclrds_max) idx_bcext_nclrds = idx_bc_nclrds_max

               ! retrieve absorption enhancement factor for within-ice BC
               enh_fct = bcenh(bnd_idx,idx_bcint_nclrds,idx_bcint_icerds)

               ! get BC optical properties (moved from above)
               ! aerosol species 1 optical properties (within-ice BC)
               ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx,idx_bcint_nclrds)
               asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx,idx_bcint_nclrds)
               ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx,idx_bcint_nclrds)*enh_fct

               ! aerosol species 2 optical properties (external BC)
               ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx,idx_bcext_nclrds)
               asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx,idx_bcext_nclrds)
               ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx,idx_bcext_nclrds)

#else
               ! bulk aerosol treatment (BC optical properties independent
               ! of BC and ice grain size)
               ! aerosol species 1 optical properties (within-ice BC)
               ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx)
               asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx)
               ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx)

               ! aerosol species 2 optical properties
               ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx)
               asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx)
               ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx)
#endif

               ! Calculate single-scattering albedo for internal mixing of dust-snow
               IF (use_dust_snow_internal_mixing) THEN
                   IF (bnd_idx < 4) THEN
                      C_dust_total = mss_cnc_aer_lcl(i,5) + mss_cnc_aer_lcl(i,6) &
                                   + mss_cnc_aer_lcl(i,7) + mss_cnc_aer_lcl(i,8)
                      C_dust_total = C_dust_total * 1.0E+06_r8
                      IF(C_dust_total > 0._r8) THEN
                         IF (flg_slr_in == 1) THEN
                            R_1_omega_tmp = dust_clear_d0(bnd_idx) &
                                          + dust_clear_d2(bnd_idx)*(C_dust_total**dust_clear_d1(bnd_idx))
                         ELSE
                            R_1_omega_tmp = dust_cloudy_d0(bnd_idx) &
                                          + dust_cloudy_d2(bnd_idx)*(C_dust_total**dust_cloudy_d1(bnd_idx))
                         ENDIF
                         ss_alb_snw_lcl(i) = 1.0_r8 - (1.0_r8 - ss_alb_snw_lcl(i)) *R_1_omega_tmp
                      ENDIF
                   ENDIF
                   DO j = 5,8,1
                      ss_alb_aer_lcl(j)        = 0._r8
                      asm_prm_aer_lcl(j)       = 0._r8
                      ext_cff_mss_aer_lcl(j)   = 0._r8
                   ENDDO
               ENDIF

               !mgf--

               L_snw(i)   = h2osno_ice_lcl(i)+h2osno_liq_lcl(i)
               tau_snw(i) = L_snw(i)*ext_cff_mss_snw_lcl(i)

               DO j=1,sno_nbr_aer
                  IF (use_dust_snow_internal_mixing .and. (j >= 5)) THEN
                    L_aer(i,j)  = 0._r8
                  ELSE
                    L_aer(i,j)   = L_snw(i)*mss_cnc_aer_lcl(i,j)
                  ENDIF
                  tau_aer(i,j) = L_aer(i,j)*ext_cff_mss_aer_lcl(j)
               ENDDO

               tau_sum   = 0._r8
               omega_sum = 0._r8
               g_sum     = 0._r8

               DO j=1,sno_nbr_aer
                  tau_sum    = tau_sum + tau_aer(i,j)
                  omega_sum  = omega_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j))
                  g_sum      = g_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j)*asm_prm_aer_lcl(j))
               ENDDO

               tau(i)    = tau_sum + tau_snw(i)
               omega(i)  = (1/tau(i))*(omega_sum+(ss_alb_snw_lcl(i)*tau_snw(i)))
               g(i)      = (1/(tau(i)*omega(i)))*(g_sum+ (asm_prm_snw_lcl(i)*ss_alb_snw_lcl(i)*tau_snw(i)))
            ENDDO ! endWeighted Mie parameters of each layer

            ! DELTA transformations, IF requested
            IF (DELTA == 1) THEN
               DO i=snl_top,snl_btm,1
                  g_star(i)     = g(i)/(1+g(i))
                  omega_star(i) = ((1-(g(i)**2))*omega(i)) / (1-(omega(i)*(g(i)**2)))
                  tau_star(i)   = (1-(omega(i)*(g(i)**2)))*tau(i)
               ENDDO
            ELSE
               DO i=snl_top,snl_btm,1
                  g_star(i)     = g(i)
                  omega_star(i) = omega(i)
                  tau_star(i)   = tau(i)
               ENDDO
            ENDIF

            ! Begin radiative transfer solver
            ! Given input vertical profiles of optical properties, evaluate the
            ! monochromatic Delta-Eddington adding-doubling solution

            ! note that trndir, trntdr, trndif, rupdir, rupdif, rdndif
            ! are variables at the layer interface,
            ! for snow with layers rangeing from snl_top to snl_btm
            ! there are snl_top to snl_btm+1 layer interface
            snl_btm_itf = snl_btm + 1

            DO i = snl_top,snl_btm_itf,1
               trndir(i) = c0
               trntdr(i) = c0
               trndif(i) = c0
               rupdir(i) = c0
               rupdif(i) = c0
               rdndif(i) = c0
            ENDDO

            ! initialize top interface of top layer
            trndir(snl_top) = c1
            trntdr(snl_top) = c1
            trndif(snl_top) = c1
            rdndif(snl_top) = c0

            ! begin main level loop
            ! for layer interfaces except for the very bottom
            DO i = snl_top,snl_btm,1

               ! initialize all layer apparent optical properties to 0
               rdir  (i) = c0
               rdif_a(i) = c0
               rdif_b(i) = c0
               tdir  (i) = c0
               tdif_a(i) = c0
               tdif_b(i) = c0
               trnlay(i) = c0

               ! compute next layer Delta-eddington solution only IF total transmission
               ! of radiation to the interface just above the layer exceeds trmin.

               IF (trntdr(i) > trmin ) THEN

                  ! calculation over layers with penetrating radiation

                  ! delta-transformed single-scattering properties
                  ! of this layer
                  ts = tau_star(i)
                  ws = omega_star(i)
                  gs = g_star(i)

                 ! Delta-Eddington solution expressions
                 ! n(uu,et)         = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et)
                 ! u(w,gg,e)        = c1p5*(c1 - w*gg)/e
                 ! el(w,gg)         = sqrt(c3*(c1-w)*(c1 - w*gg))
                 lm = sqrt(c3*(c1-ws)*(c1 - ws*gs))  !lm = el(ws,gs)
                 ue = c1p5*(c1 - ws*gs)/lm           !ue = u(ws,gs,lm)
                 extins = max(exp_min, exp(-lm*ts))
                 ne = ((ue+c1)*(ue+c1)/extins) - ((ue-c1)*(ue-c1)*extins) !ne = n(ue,extins)

                 ! first calculation of rdif, tdif using Delta-Eddington formulas
                 ! rdif_a(k) = (ue+c1)*(ue-c1)*(c1/extins - extins)/ne
                 rdif_a(i) = (ue**2-c1)*(c1/extins - extins)/ne
                 tdif_a(i) = c4*ue/ne

                 ! evaluate rdir,tdir for direct beam
                 trnlay(i) = max(exp_min, exp(-ts/mu_not))

                 ! Delta-Eddington solution expressions
                 ! alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
                 ! agamm(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
                 ! alp = alpha(ws,mu_not,gs,lm)
                 ! gam = agamm(ws,mu_not,gs,lm)
                 alp = cp75*ws*mu_not*((c1 + gs*(c1-ws))/(c1 - lm*lm*mu_not*mu_not))
                 gam = cp5*ws*((c1 + c3*gs*(c1-ws)*mu_not*mu_not)/(c1-lm*lm*mu_not*mu_not))
                 apg = alp + gam
                 amg = alp - gam

                 rdir(i) = apg*rdif_a(i) +  amg*(tdif_a(i)*trnlay(i) - c1)
                 tdir(i) = apg*tdif_a(i) + (amg* rdif_a(i)-apg+c1)*trnlay(i)

                 ! recalculate rdif,tdif using direct angular integration over rdir,tdir,
                 ! since Delta-Eddington rdif formula is not well-behaved (it is usually
                 ! biased low and can even be negative); USE ngmax angles and gaussian
                 ! integration for most accuracy:
                 R1 = rdif_a(i) ! USE R1 as temporary
                 T1 = tdif_a(i) ! USE T1 as temporary
                 swt = c0
                 smr = c0
                 smt = c0
                 DO ng=1,ngmax
                    mu  = difgauspt(ng)
                    gwt = difgauswt(ng)
                    swt = swt + mu*gwt
                    trn = max(exp_min, exp(-ts/mu))
                    ! alp = alpha(ws,mu,gs,lm)
                    ! gam = agamm(ws,mu,gs,lm)
                    alp = cp75*ws*mu*((c1 + gs*(c1-ws))/(c1 - lm*lm*mu*mu))
                    gam = cp5*ws*((c1 + c3*gs*(c1-ws)*mu*mu)/(c1-lm*lm*mu*mu))
                    apg = alp + gam
                    amg = alp - gam
                    rdr = apg*R1 + amg*T1*trn - amg
                    tdr = apg*T1 + amg*R1*trn - apg*trn + trn
                    smr = smr + mu*rdr*gwt
                    smt = smt + mu*tdr*gwt
                 ENDDO      ! ng
                 rdif_a(i) = smr/swt
                 tdif_a(i) = smt/swt

                 ! homogeneous layer
                 rdif_b(i) = rdif_a(i)
                 tdif_b(i) = tdif_a(i)

               ENDIF ! trntdr(k) > trmin

               ! Calculate the solar beam transmission, total transmission, and
               ! reflectivity for diffuse radiation from below at interface i,
               ! the top of the current layer k:
               !
               !              layers       interface
               !
               !       ---------------------  i-1
               !                i-1
               !       ---------------------  i
               !                 i
               !       ---------------------

               trndir(i+1) = trndir(i)*trnlay(i)
               refkm1      = c1/(c1 - rdndif(i)*rdif_a(i))
               tdrrdir     = trndir(i)*rdir(i)
               tdndif      = trntdr(i) - trndir(i)
               trntdr(i+1) = trndir(i)*tdir(i) + &
                    (tdndif + tdrrdir*rdndif(i))*refkm1*tdif_a(i)
               rdndif(i+1) = rdif_b(i) + &
                    (tdif_b(i)*rdndif(i)*refkm1*tdif_a(i))
               trndif(i+1) = trndif(i)*refkm1*tdif_a(i)

            ENDDO       !   END main level loop


            ! compute reflectivity to direct and diffuse radiation for layers
            ! below by adding succesive layers starting from the underlying
            ! ground and working upwards:
            !
            !              layers       interface
            !
            !       ---------------------  i
            !                 i
            !       ---------------------  i+1
            !                i+1
            !       ---------------------

            ! set the underlying ground albedo == albedo of near-IR
            ! unless bnd_idx == 1, for visible
            rupdir(snl_btm_itf) = albsfc(2)
            rupdif(snl_btm_itf) = albsfc(2)
            IF (bnd_idx == 1) THEN
                rupdir(snl_btm_itf) = albsfc(1)
                rupdif(snl_btm_itf) = albsfc(1)
            ENDIF

            DO i=snl_btm,snl_top,-1
               ! interface scattering
               refkp1        = c1/( c1 - rdif_b(i)*rupdif(i+1))
               ! dir from top layer plus exp tran ref from lower layer, interface
               ! scattered and tran thru top layer from below, plus diff tran ref
               ! from lower layer with interface scattering tran thru top from below
               rupdir(i) = rdir(i) &
                    + (        trnlay(i)  *rupdir(i+1) &
                    +  (tdir(i)-trnlay(i))*rupdif(i+1))*refkp1*tdif_b(i)
               ! dif from top layer from above, plus dif tran upwards reflected and
               ! interface scattered which tran top from below
               rupdif(i) = rdif_a(i) + tdif_a(i)*rupdif(i+1)*refkp1*tdif_b(i)
            ENDDO       ! i

            ! net flux (down-up) at each layer interface from the
            ! snow top (i = snl_top) to bottom interface above land (i = snl_btm_itf)
            ! the interface reflectivities and transmissivities required
            ! to evaluate interface fluxes are returned from solution_dEdd;
            ! now compute up and down fluxes for each interface, using the
            ! combined layer properties at each interface:
            !
            !              layers       interface
            !
            !       ---------------------  i
            !                 i
            !       ---------------------

            DO i = snl_top, snl_btm_itf
               ! interface scattering
               refk          = c1/(c1 - rdndif(i)*rupdif(i))
               ! dir tran ref from below times interface scattering, plus diff
               ! tran and ref from below times interface scattering
               ! fdirup(i) = (trndir(i)*rupdir(i) + &
               !                 (trntdr(i)-trndir(i))  &
               !                 *rupdif(i))*refk
               ! dir tran plus total diff trans times interface scattering plus
               ! dir tran with up dir ref and down dif ref times interface scattering
               ! fdirdn(i) = trndir(i) + (trntdr(i) &
               !               - trndir(i) + trndir(i)  &
               !               *rupdir(i)*rdndif(i))*refk
               ! diffuse tran ref from below times interface scattering
               ! fdifup(i) = trndif(i)*rupdif(i)*refk
               ! diffuse tran times interface scattering
               ! fdifdn(i) = trndif(i)*refk

               ! netflux, down - up
               ! dfdir = fdirdn - fdirup
               dfdir(i) = trndir(i) &
                           + (trntdr(i)-trndir(i)) * (c1 - rupdif(i)) * refk &
                           -  trndir(i)*rupdir(i)  * (c1 - rdndif(i)) * refk
               IF (dfdir(i) < puny) dfdir(i) = c0
               ! dfdif = fdifdn - fdifup
               dfdif(i) = trndif(i) * (c1 - rupdif(i)) * refk
               IF (dfdif(i) < puny) dfdif(i) = c0
            ENDDO       ! k

            ! SNICAR_AD_RT is called twice for direct and diffuse incident fluxes
            ! direct incident
            IF (flg_slr_in == 1) THEN
              albedo = rupdir(snl_top)
              dftmp  = dfdir
              refk   = c1/(c1 - rdndif(snl_top)*rupdif(snl_top))
              F_sfc_pls = (trndir(snl_top)*rupdir(snl_top) + &
                          (trntdr(snl_top)-trndir(snl_top))  &
                           *rupdif(snl_top))*refk
            !diffuse incident
            ELSE
              albedo = rupdif(snl_top)
              dftmp  = dfdif
              refk   = c1/(c1 - rdndif(snl_top)*rupdif(snl_top))
              F_sfc_pls = trndif(snl_top)*rupdif(snl_top)*refk
            ENDIF

            ! Absorbed flux in each layer
            DO i=snl_top,snl_btm,1
              F_abs(i) = dftmp(i)-dftmp(i+1)
              flx_abs_lcl(i,bnd_idx) = F_abs(i)

              ! ERROR check: negative absorption
              IF (p_is_master) THEN
                 IF (flx_abs_lcl(i,bnd_idx) < -0.00001) THEN
                    write (iulog,"(a,e13.6,a,i6)") "SNICAR ERROR: negative absoption : ", flx_abs_lcl(i,bnd_idx)
                    write(iulog,*) "SNICAR_AD STATS: snw_rds(0)= ", snw_rds(0)
                    write(iulog,*) "SNICAR_AD STATS: L_snw(0)= ", L_snw(0)
                    write(iulog,*) "SNICAR_AD STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
                    write(iulog,*) "SNICAR_AD STATS: soot1(0)= ", mss_cnc_aer_lcl(0,1)
                    write(iulog,*) "SNICAR_AD STATS: soot2(0)= ", mss_cnc_aer_lcl(0,2)
                    write(iulog,*) "SNICAR_AD STATS: dust1(0)= ", mss_cnc_aer_lcl(0,3)
                    write(iulog,*) "SNICAR_AD STATS: dust2(0)= ", mss_cnc_aer_lcl(0,4)
                    write(iulog,*) "SNICAR_AD STATS: dust3(0)= ", mss_cnc_aer_lcl(0,5)
                    write(iulog,*) "SNICAR_AD STATS: dust4(0)= ", mss_cnc_aer_lcl(0,6)
                    CALL abort
                 ENDIF
              ENDIF
            ENDDO

            ! absobed flux by the underlying ground
            F_btm_net = dftmp(snl_btm_itf)

            ! note here, snl_btm_itf = 1 by snow column set up in CLM
            flx_abs_lcl(1,bnd_idx) = F_btm_net

            IF (flg_nosnl == 1) THEN
               ! If there are no snow layers (but still snow), all absorbed energy must be in top soil layer
               !flx_abs_lcl(:,bnd_idx) = 0._r8
               !flx_abs_lcl(1,bnd_idx) = F_abs(0) + F_btm_net

               ! changed on 20070408:
               ! OK to put absorbed energy in the fictitous snow layer because routine SurfaceRadiation
               ! handles the CASE of no snow layers. Then, IF a snow layer is addded between now and
               ! SurfaceRadiation (called in CanopyHydrology), absorbed energy will be properly distributed.
               flx_abs_lcl(0,bnd_idx) = F_abs(0)
               flx_abs_lcl(1,bnd_idx) = F_btm_net
            ENDIF

            !Underflow check (we've already tripped the error condition above)
            DO i=snl_top,1,1
               IF (flx_abs_lcl(i,bnd_idx) < 0._r8) THEN
                  flx_abs_lcl(i,bnd_idx) = 0._r8
               ENDIF
            ENDDO

            F_abs_sum = 0._r8
            DO i=snl_top,snl_btm,1
               F_abs_sum = F_abs_sum + F_abs(i)
            ENDDO

            !ENDDO !ENDDO WHILE (flg_dover > 0)

            ! Energy conservation check:
            ! Incident direct+diffuse radiation equals (absorbed+bulk_transmitted+bulk_reflected)
            energy_sum = (mu_not*pi*flx_slrd_lcl(bnd_idx)) + flx_slri_lcl(bnd_idx) - (F_abs_sum + F_btm_net + F_sfc_pls)
            IF (p_is_master) THEN
               IF (abs(energy_sum) > 0.00001_r8) THEN
                  write (iulog,"(a,e13.6,a,i6)") "SNICAR ERROR: Energy conservation error of : ", energy_sum
                  write(iulog,*) "F_abs_sum: ",F_abs_sum
                  write(iulog,*) "F_btm_net: ",F_btm_net
                  write(iulog,*) "F_sfc_pls: ",F_sfc_pls
                  write(iulog,*) "mu_not*pi*flx_slrd_lcl(bnd_idx): ", mu_not*pi*flx_slrd_lcl(bnd_idx)
                  write(iulog,*) "flx_slri_lcl(bnd_idx)", flx_slri_lcl(bnd_idx)
                  write(iulog,*) "bnd_idx", bnd_idx
                  write(iulog,*) "F_abs", F_abs
                  write(iulog,*) "albedo", albedo
                  CALL abort
               ENDIF
            ENDIF

            albout_lcl(bnd_idx) = albedo
            ! Check that albedo is less than 1
            IF (p_is_master) THEN
               IF (albout_lcl(bnd_idx) > 1.0) THEN
                  write (iulog,*) "SNICAR ERROR: Albedo > 1.0: "
                  write (iulog,*) "SNICAR STATS: bnd_idx= ",bnd_idx
                  write (iulog,*) "SNICAR STATS: albout_lcl(bnd)= ",albout_lcl(bnd_idx), &
                     " albsfc_lcl(bnd_idx)= ",albsfc_lcl(bnd_idx)
                  write (iulog,*) "SNICAR STATS: h2osno= ", h2osno_lcl, " snl= ", snl_lcl
                  write (iulog,*) "SNICAR STATS: coszen= ", coszen, " flg_slr= ", flg_slr_in

                  write (iulog,*) "SNICAR STATS: soot(-4)= ", mss_cnc_aer_lcl(-4,1)
                  write (iulog,*) "SNICAR STATS: soot(-3)= ", mss_cnc_aer_lcl(-3,1)
                  write (iulog,*) "SNICAR STATS: soot(-2)= ", mss_cnc_aer_lcl(-2,1)
                  write (iulog,*) "SNICAR STATS: soot(-1)= ", mss_cnc_aer_lcl(-1,1)
                  write (iulog,*) "SNICAR STATS: soot(0)= ", mss_cnc_aer_lcl(0,1)

                  write (iulog,*) "SNICAR STATS: L_snw(-4)= ", L_snw(-4)
                  write (iulog,*) "SNICAR STATS: L_snw(-3)= ", L_snw(-3)
                  write (iulog,*) "SNICAR STATS: L_snw(-2)= ", L_snw(-2)
                  write (iulog,*) "SNICAR STATS: L_snw(-1)= ", L_snw(-1)
                  write (iulog,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)

                  write (iulog,*) "SNICAR STATS: snw_rds(-4)= ", snw_rds(-4)
                  write (iulog,*) "SNICAR STATS: snw_rds(-3)= ", snw_rds(-3)
                  write (iulog,*) "SNICAR STATS: snw_rds(-2)= ", snw_rds(-2)
                  write (iulog,*) "SNICAR STATS: snw_rds(-1)= ", snw_rds(-1)
                  write (iulog,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(0)

                  CALL abort
               ENDIF
            ENDIF

         ENDDO   ! loop over wvl bands


         ! Weight output NIR albedo appropriately
         albout(1) = albout_lcl(1)
         flx_sum         = 0._r8
         DO bnd_idx= nir_bnd_bgn,nir_bnd_end
            flx_sum = flx_sum + flx_wgt(bnd_idx)*albout_lcl(bnd_idx)
         ENDDO
         albout(2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))

         ! Weight output NIR absorbed layer fluxes (flx_abs) appropriately
         flx_abs(:,1) = flx_abs_lcl(:,1)
         DO i=snl_top,1,1
            flx_sum = 0._r8
            DO bnd_idx= nir_bnd_bgn,nir_bnd_end
               flx_sum = flx_sum + flx_wgt(bnd_idx)*flx_abs_lcl(i,bnd_idx)
            ENDDO
            flx_abs(i,2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
         ENDDO

         ! near-IR direct albedo/absorption adjustment for high solar zenith angles
         ! solar zenith angle parameterization
         ! calculate the scaling factor for NIR direct albedo IF SZA>75 degree
         IF ((mu_not < mu_75) .and. (flg_slr_in == 1)) THEN
            sza_c1 = sza_a0 + sza_a1 * mu_not + sza_a2 * mu_not**2
            sza_c0 = sza_b0 + sza_b1 * mu_not + sza_b2 * mu_not**2
            sza_factor = sza_c1 * (log10(snw_rds_lcl(snl_top) * c1) - c6) + sza_c0
            flx_sza_adjust  = albout(2) * (sza_factor-c1) * sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
            albout(2) = albout(2) * sza_factor
            flx_abs(snl_top,2) = flx_abs(snl_top,2) - flx_sza_adjust
         ENDIF

      ! If snow < minimum_snow, but > 0, and there is sun, set albedo to underlying surface albedo
      ELSEIF ( (coszen > 0._r8) .and. (h2osno_lcl < min_snw) .and. (h2osno_lcl > 0._r8) ) THEN
         albout(1) = albsfc(1)
         albout(2) = albsfc(2)

         ! There is either zero snow, or no sun
      ELSE
         albout(1) = 0._r8
         albout(2) = 0._r8
      ENDIF    ! IF column has snow and coszen > 0

   ! END associate

   END SUBROUTINE SNICAR_AD_RT
   !-----------------------------------------------------------------------


   SUBROUTINE SnowAge_grain( dtime          , snl            , dz              ,&
                             qflx_snow_grnd , qflx_snwcp_ice , qflx_snofrz_lyr ,&
                             do_capsnow     , frac_sno       , h2osno          ,&
                             h2osno_liq     , h2osno_ice     , t_soisno        ,&
                             t_grnd         , forc_t         , snw_rds          )
!-----------------------------------------------------------------------
!
! !DESCRIPTION:
!  Updates the snow effective grain size (radius).
!  Contributions to grain size evolution are from:
!    1. vapor redistribution (dry snow)
!    2. liquid water redistribution (wet snow)
!    3. re-freezing of liquid water
!
!  Vapor redistribution: Method is to retrieve 3 best-bit parameters that
!  depend on snow temperature, temperature gradient, and density,
!  that are derived from the microphysical model described in:
!  Flanner and Zender (2006), Linking snowpack microphysics and albedo
!  evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834.
!  The parametric equation has the form:
!  dr/dt = drdt_0*(tau/(dr_fresh+tau))^(1/kappa), WHERE:
!    r is the effective radius,
!    tau and kappa are best-fit parameters,
!    drdt_0 is the initial rate of change of effective radius, and
!    dr_fresh is the difference between the current and fresh snow states
!   (r_current - r_fresh).
!
!  Liquid water redistribution: Apply the grain growth FUNCTION from:
!    Brun, E. (1989), Investigation of wet-snow metamorphism in respect of
!    liquid-water content, Annals of Glaciology, 13, 22-26.
!    There are two parameters that describe the grain growth rate as
!    a FUNCTION of snow liquid water content (LWC). The "LWC=0" parameter
!    is zeroed here because we are accounting for dry snowing with a
!    different representation
!
!  Re-freezing of liquid water: Assume that re-frozen liquid water clumps
!    into an arbitrarily large effective grain size (snw_rds_refrz).
!    The phenomenon is observed (Grenfell), but so far unquantified, as far as
!    I am aware.
!
! !USES:
!
! DAI, Dec. 29, 2022
!-----------------------------------------------------------------------
! !ARGUMENTS:

   IMPLICIT NONE

   real(r8) , intent(in)   ::  dtime                           ! land model time step [sec]

   integer  , intent(in)   ::  snl                             ! negative number of snow layers (col) [nbr]
   real(r8) , intent(in)   ::  dz               ( maxsnl+1:1 ) ! layer thickness (col,lyr) [m]

   real(r8) , intent(in)   ::  qflx_snow_grnd                  ! snow on ground after interception (col) [kg m-2 s-1]
   real(r8) , intent(in)   ::  qflx_snwcp_ice                  ! excess precipitation due to snow capping [kg m-2 s-1]
   real(r8) , intent(in)   ::  qflx_snofrz_lyr  ( maxsnl+1:0 ) ! snow freezing rate (col,lyr) [kg m-2 s-1]

   logical  , intent(in)   ::  do_capsnow                      ! true => DO snow capping
   real(r8) , intent(in)   ::  frac_sno                        ! fraction of ground covered by snow (0 to 1)
   real(r8) , intent(in)   ::  h2osno                          ! snow water (col) [mm H2O]
   real(r8) , intent(in)   ::  h2osno_liq       ( maxsnl+1:0 ) ! liquid water content (col,lyr) [kg m-2]
   real(r8) , intent(in)   ::  h2osno_ice       ( maxsnl+1:0 ) ! ice content (col,lyr) [kg m-2]

   real(r8) , intent(in)   ::  t_soisno         ( maxsnl+1:1 ) ! soil and snow temperature (col,lyr) [K]
   real(r8) , intent(in)   ::  t_grnd                          ! ground temperature (col) [K]
   real(r8) , intent(in)   ::  forc_t                          ! Atmospheric temperature (col) [K]

   real(r8) , intent(inout)  ::  snw_rds        ( maxsnl+1:0 ) ! effective grain radius (col,lyr) [microns, m-6]

   ! DAI, Dec. 29, 2022
   !-----------------------------------------------------------------------
   !
   ! !LOCAL VARIABLES:
   integer :: snl_top          ! top snow layer index [idx]
   integer :: snl_btm          ! bottom snow layer index [idx]
   integer :: i                ! layer index [idx]
                               ! integer :: c_idx                                                 ! column index [idx]
   integer :: fc               ! snow column filter index [idx]
   integer :: T_idx            ! snow aging lookup table temperature index [idx]
   integer :: Tgrd_idx         ! snow aging lookup table temperature gradient index [idx]
   integer :: rhos_idx         ! snow aging lookup table snow density index [idx]
   real(r8) :: t_snotop        ! temperature at upper layer boundary [K]
   real(r8) :: t_snobtm        ! temperature at lower layer boundary [K]
   real(r8) :: dTdz(maxsnl:0)  ! snow temperature gradient (col,lyr) [K m-1]
   real(r8) :: bst_tau         ! snow aging parameter retrieved from lookup table [hour]
   real(r8) :: bst_kappa       ! snow aging parameter retrieved from lookup table [unitless]
   real(r8) :: bst_drdt0       ! snow aging parameter retrieved from lookup table [um hr-1]
   real(r8) :: dr              ! incremental change in snow effective radius [um]
   real(r8) :: dr_wet          ! incremental change in snow effective radius from wet growth [um]
   real(r8) :: dr_fresh        ! difference between fresh snow r_e and current r_e [um]
   real(r8) :: newsnow         ! fresh snowfall [kg m-2]
   real(r8) :: refrzsnow       ! re-frozen snow [kg m-2]
   real(r8) :: frc_newsnow     ! fraction of layer mass that is new snow [frc]
   real(r8) :: frc_oldsnow     ! fraction of layer mass that is old snow [frc]
   real(r8) :: frc_refrz       ! fraction of layer mass that is re-frozen snow [frc]
   real(r8) :: frc_liq         ! fraction of layer mass that is liquid water[frc]
   real(r8) :: rhos            ! snow density [kg m-3]
   real(r8) :: h2osno_lyr      ! liquid + solid H2O in snow layer [kg m-2]
   real(r8) :: cdz(maxsnl+1:0) ! column average layer thickness [m]
   real(r8) :: snw_rds_fresh   ! fresh snow radius [microns]

   real(r8) :: snot_top        ! temperature in top snow layer (col) [K]
   real(r8) :: dTdz_top        ! temperature gradient in top layer (col) [K m-1]
   real(r8) :: snw_rds_top     ! effective grain radius, top layer (col) [microns, m-6]
   real(r8) :: sno_liq_top     ! liquid water fraction (mass) in top snow layer (col) [frc]

   !--------------------------------------------------------------------------!

 !  associate(                                          &
 !       snl             => col_pp%snl             , & ! Input:  [integer  (:)   ]  negative number of snow layers (col) [nbr]
 !       dz              => col_pp%dz              , & ! Input:  [real(r8) (:,:) ]  layer thickness (col,lyr) [m]

 !       qflx_snow_grnd  => col_wf%qflx_snow_grnd  , & ! Input:  [real(r8) (:)   ]  snow on ground after interception (col) [kg m-2 s-1]
 !       qflx_snwcp_ice  => col_wf%qflx_snwcp_ice  , & ! Input:  [real(r8) (:)   ]  excess precipitation due to snow capping [kg m-2 s-1]
 !       qflx_snofrz_lyr => col_wf%qflx_snofrz_lyr , & ! Input:  [real(r8) (:,:) ]  snow freezing rate (col,lyr) [kg m-2 s-1]

 !       do_capsnow      => col_ws%do_capsnow      , & ! Input:  [logical  (:)   ]  true => DO snow capping
 !       frac_sno        => col_ws%frac_sno_eff    , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by snow (0 to 1)
 !       h2osno          => col_ws%h2osno          , & ! Input:  [real(r8) (:)   ]  snow water (col) [mm H2O]
 !       h2osno_liq      => col_ws%h2osno_liq      , & ! Input:  [real(r8) (:,:) ]  liquid water content (col,lyr) [kg m-2]
 !       h2osno_ice      => col_ws%h2osno_ice      , & ! Input:  [real(r8) (:,:) ]  ice content (col,lyr) [kg m-2]
 !       snw_rds         => col_ws%snw_rds         , & ! Output: [real(r8) (:,:) ]  effective grain radius (col,lyr) [microns, m-6]
 !       snw_rds_top     => col_ws%snw_rds_top     , & ! Output: [real(r8) (:)   ]  effective grain radius, top layer (col) [microns, m-6]
 !       sno_liq_top     => col_ws%sno_liq_top     , & ! Output: [real(r8) (:)   ]  liquid water fraction (mass) in top snow layer (col) [frc]

 !       t_soisno        => col_es%t_soisno        , & ! Input:  [real(r8) (:,:) ]  soil and snow temperature (col,lyr) [K]
 !       t_grnd          => col_es%t_grnd          , & ! Input:  [real(r8) (:)   ]  ground temperature (col) [K]
 !       snot_top        => col_es%snot_top        , & ! Output: [real(r8) (:)   ]  temperature in top snow layer (col) [K]
 !       dTdz_top        => col_es%dTdz_top          & ! Output: [real(r8) (:)   ]  temperature gradient in top layer (col) [K m-1]
 !       )


      IF (snl < 0 .and. h2osno > 0._r8) THEN

         snl_btm = 0
         snl_top = snl + 1

         cdz(snl_top:snl_btm)=frac_sno*dz(snl_top:snl_btm)

         ! loop over snow layers
         DO i = snl_top, snl_btm, 1
            !
            !**********  1. DRY SNOW AGING  ***********
            !
            h2osno_lyr = h2osno_liq(i) + h2osno_ice(i)

            ! temperature gradient
            IF (i == snl_top) THEN
               ! top layer
               t_snotop = t_soisno(snl_top)
               t_snobtm = (t_soisno(i+1)*dz(i) &
                         + t_soisno(i)*dz(i+1)) &
                         / (dz(i)+dz(i+1))
            ELSE
               t_snotop = (t_soisno(i-1)*dz(i) &
                         + t_soisno(i)*dz(i-1)) &
                         / (dz(i)+dz(i-1))
               t_snobtm = (t_soisno(i+1)*dz(i) &
                         + t_soisno(i)*dz(i+1)) &
                         / (dz(i)+dz(i+1))
            ENDIF

            dTdz(i) = abs((t_snotop - t_snobtm) / cdz(i))

            ! snow density
            rhos = (h2osno_liq(i)+h2osno_ice(i)) / cdz(i)

            ! make sure rhos doesn't drop below 50 (see rhos_idx below)
            rhos=max(50._r8,rhos)

            ! best-fit table indecies
            T_idx    = nint((t_soisno(i)-223) / 5) + 1
            Tgrd_idx = nint(dTdz(i) / 10) + 1
            rhos_idx = nint((rhos-50) / 50) + 1

            ! boundary check:
            IF (T_idx < idx_T_min) THEN
               T_idx = idx_T_min
            ENDIF
            IF (T_idx > idx_T_max) THEN
               T_idx = idx_T_max
            ENDIF
            IF (Tgrd_idx < idx_Tgrd_min) THEN
               Tgrd_idx = idx_Tgrd_min
            ENDIF
            IF (Tgrd_idx > idx_Tgrd_max) THEN
               Tgrd_idx = idx_Tgrd_max
            ENDIF
            IF (rhos_idx < idx_rhos_min) THEN
               rhos_idx = idx_rhos_min
            ENDIF
            IF (rhos_idx > idx_rhos_max) THEN
               rhos_idx = idx_rhos_max
            ENDIF

            ! best-fit parameters
            bst_tau   = snowage_tau(rhos_idx,Tgrd_idx,T_idx)
            bst_kappa = snowage_kappa(rhos_idx,Tgrd_idx,T_idx)
            bst_drdt0 = snowage_drdt0(rhos_idx,Tgrd_idx,T_idx)

            ! change in snow effective radius, using best-fit parameters
            ! added checks suggested by mgf. --HW 10/15/2015
            dr_fresh = snw_rds(i)-snw_rds_min

#ifdef MODAL_AER
            IF ( abs(dr_fresh) < 1.0e-8_r8 ) THEN
               dr_fresh = 0.0_r8
            ELSEIF ( dr_fresh < 0.0_r8 ) THEN
               IF (p_is_master) THEN
                  write(iulog,*) "dr_fresh = ", dr_fresh, snw_rds(i), snw_rds_min
                  CALL abort
               ENDIF
            ENDIF

            dr = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))**(1._r8/bst_kappa)) * (dtime/3600._r8)
#else
            dr = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))**(1/bst_kappa)) * (dtime/3600)
#endif
            !
            !**********  2. WET SNOW AGING  ***********
            !
            ! We are assuming wet and dry evolution occur simultaneously, and
            ! the contributions from both can be summed.
            ! This is justified by setting the linear offset constant C1_liq_Brun89 to zero [Brun, 1989]

            ! liquid water faction
            frc_liq = min(0.1_r8, (h2osno_liq(i) / (h2osno_liq(i)+h2osno_ice(i))))

            !dr_wet = 1E6_r8*(dtime*(C1_liq_Brun89 + C2_liq_Brun89*(frc_liq**(3))) / (4*SHR_CONST_PI*(snw_rds(i)/1E6)**(2)))
            !simplified, units of microns:
            dr_wet = 1E18_r8*(dtime*(C2_liq_Brun89*(frc_liq**(3))) / (4*SHR_CONST_PI*snw_rds(i)**(2)))

            dr = dr + dr_wet

            !
            !**********  3. SNOWAGE SCALING (TURNED OFF BY DEFAULT)  *************
            !
            ! Multiply rate of change of effective radius by some constant, xdrdt
            IF (flg_snoage_scl) THEN
               dr = dr*xdrdt
            ENDIF

            !
            !**********  4. INCREMENT EFFECTIVE RADIUS, ACCOUNTING FOR:  ***********
            !               DRY AGING
            !               WET AGING
            !               FRESH SNOW
            !               RE-FREEZING
            !
            ! new snowfall [kg/m2]
            IF (do_capsnow .and. .not. use_extrasnowlayers) THEN
               newsnow = max(0._r8, (qflx_snwcp_ice*dtime))
            ELSE
               newsnow = max(0._r8, (qflx_snow_grnd*dtime))
            ENDIF

            ! snow that has re-frozen [kg/m2]
            refrzsnow = max(0._r8, (qflx_snofrz_lyr(i)*dtime))

            ! fraction of layer mass that is re-frozen
            frc_refrz = refrzsnow / h2osno_lyr

            ! fraction of layer mass that is new snow
            IF (i == snl_top) THEN
               frc_newsnow = newsnow / h2osno_lyr
            ELSE
               frc_newsnow = 0._r8
            ENDIF

            IF ((frc_refrz + frc_newsnow) > 1._r8) THEN
               frc_refrz = frc_refrz / (frc_refrz + frc_newsnow)
               frc_newsnow = 1._r8 - frc_refrz
               frc_oldsnow = 0._r8
            ELSE
               frc_oldsnow = 1._r8 - frc_refrz - frc_newsnow
            ENDIF

            ! temperature dependent fresh grain size
            snw_rds_fresh = FreshSnowRadius (forc_t)

            ! mass-weighted mean of fresh snow, old snow, and re-frozen snow effective radius
            snw_rds(i) = (snw_rds(i)+dr)*frc_oldsnow + snw_rds_fresh*frc_newsnow + snw_rds_refrz*frc_refrz
            !
            !**********  5. CHECK BOUNDARIES   ***********
            !
            ! boundary check
            IF (snw_rds(i) < snw_rds_min) THEN
               snw_rds(i) = snw_rds_min
            ENDIF

            IF (snw_rds(i) > snw_rds_max) THEN
               snw_rds(i) = snw_rds_max
            ENDIF

            ! set top layer variables for history files
            IF (i == snl_top) THEN
               snot_top    = t_soisno(i)
               dTdz_top    = dTdz(i)
               snw_rds_top = snw_rds(i)
               sno_liq_top = h2osno_liq(i) / (h2osno_liq(i)+h2osno_ice(i))
            ENDIF

         ENDDO
      ENDIF   ! ENDIF (snl < 0 )

      ! Special CASE: snow on ground, but not enough to have defined a snow layer:
      !   set snw_rds to fresh snow grain size:

      IF (snl >= 0 .and. h2osno > 0._r8) THEN
          snw_rds(0) = snw_rds_min
      ENDIF

   ! END associate

   END SUBROUTINE SnowAge_grain
   !-----------------------------------------------------------------------


   SUBROUTINE SnowOptics_init( fsnowoptics )

   USE MOD_NetCDFSerial

   IMPLICIT NONE

   character(len=256), intent(in) :: fsnowoptics     ! snow optical properties file name
   character(len= 32) :: subname = 'SnowOptics_init' ! SUBROUTINE name
   integer :: atm_type_index                         ! index for atmospheric type

   logical :: readvar  ! determine IF variable was read from NetCDF file
   !-----------------------------------------------------------------------

      readvar = .true.

      atm_type_index = atm_type_default
      ! Define atmospheric type
      IF (trim(snicar_atm_type) == 'default') THEN
         atm_type_index = atm_type_default
      ELSEIF (trim(snicar_atm_type) == 'mid-latitude_winter') THEN
         atm_type_index = atm_type_mid_latitude_winter
      ELSEIF (trim(snicar_atm_type) == 'mid-latitude_summer') THEN
         atm_type_index = atm_type_mid_latitude_summer
      ELSEIF (trim(snicar_atm_type) == 'sub-Arctic_winter') THEN
         atm_type_index = atm_type_sub_Arctic_winter
      ELSEIF (trim(snicar_atm_type) == 'sub-Arctic_summer') THEN
         atm_type_index = atm_type_sub_Arctic_summer
      ELSEIF (trim(snicar_atm_type) == 'summit_Greenland') THEN
         atm_type_index = atm_type_summit_Greenland
      ELSEIF (trim(snicar_atm_type) == 'high_mountain') THEN
         atm_type_index = atm_type_high_mountain
      ELSE
         IF (p_is_master) THEN
            write(iulog,*) "snicar_atm_type = ", snicar_atm_type
            CALL abort
         ENDIF
      ENDIF

      !
      ! Open optics file:
      IF (p_is_master) THEN
         write(iulog,*) 'Attempting to read snow optical properties .....'
         write(iulog,*) subname,trim(fsnowoptics)
      ENDIF

      ! direct-beam snow Mie parameters:
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_ice_drc', ss_alb_snw_drc)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_ice_drc', asm_prm_snw_drc)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_ice_drc', ext_cff_mss_snw_drc)

      ! diffuse snow Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_ice_dfs', ss_alb_snw_dfs)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_ice_dfs', asm_prm_snw_dfs)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_ice_dfs', ext_cff_mss_snw_dfs)

      !!! Direct and diffuse flux under different atmospheric conditions
      ! Direct-beam incident spectral flux:
      CALL ncio_read_bcast_serial (fsnowoptics, 'flx_wgt_dir', flx_wgt_dir)

      ! Diffuse incident spectral flux:
      CALL ncio_read_bcast_serial (fsnowoptics, 'flx_wgt_dif', flx_wgt_dif)

#ifdef MODAL_AER
      ! size-dependent BC parameters and BC enhancement factors
      IF (p_is_master) THEN
         write(iulog,*) 'Attempting to read optical properties for within-ice BC (modal aerosol treatment) ...'
      ENDIF
      !
      ! BC species 1 Mie parameters
      !
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_bc_mam', ss_alb_bc1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_bc_mam', asm_prm_bc1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_bc_mam', ext_cff_mss_bc1)
      !
      ! BC species 2 Mie parameters (identical, before enhancement factors applied)
      !
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_bc_mam', ss_alb_bc2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_bc_mam', asm_prm_bc2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_bc_mam', ext_cff_mss_bc2)
      !
      ! size-dependent BC absorption enhancement factors for within-ice BC
      CALL ncio_read_bcast_serial (fsnowoptics, 'bcint_enh_mam', bcenh)
      !
#else
      ! bulk aerosol treatment
      ! BC species 1 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_bcphil', ss_alb_bc1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_bcphil', asm_prm_bc1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_bcphil', ext_cff_mss_bc1)

      !
      ! BC species 2 Mie parameters
      !
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_bcphob', ss_alb_bc2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_bcphob', asm_prm_bc2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_bcphob', ext_cff_mss_bc2)
      !
#endif
      !
      ! OC species 1 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_ocphil', ss_alb_oc1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_ocphil', asm_prm_oc1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_ocphil', ext_cff_mss_oc1)
      !
      ! OC species 2 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_ocphob', ss_alb_oc2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_ocphob', asm_prm_oc2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_ocphob', ext_cff_mss_oc2)
      !
      ! dust species 1 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_dust01', ss_alb_dst1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_dust01', asm_prm_dst1)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_dust01', ext_cff_mss_dst1)
      !
      ! dust species 2 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_dust02', ss_alb_dst2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_dust02', asm_prm_dst2)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_dust02', ext_cff_mss_dst2)
      !
      ! dust species 3 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_dust03', ss_alb_dst3)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_dust03', asm_prm_dst3)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_dust03', ext_cff_mss_dst3)
      !
      ! dust species 4 Mie parameters
      CALL ncio_read_bcast_serial (fsnowoptics, 'ss_alb_dust04', ss_alb_dst4)
      CALL ncio_read_bcast_serial (fsnowoptics, 'asm_prm_dust04', asm_prm_dst4)
      CALL ncio_read_bcast_serial (fsnowoptics, 'ext_cff_mss_dust04', ext_cff_mss_dst4)
      !
      !

      IF (p_is_master) THEN
         write(iulog,*) 'Successfully read snow optical properties'
      ENDIF


      ! print some diagnostics:
      IF (p_is_master) THEN
         write (iulog,*) 'SNICAR: Mie single scatter albedos for direct-beam ice, rds=100um: ', &
            ss_alb_snw_drc(71,1), ss_alb_snw_drc(71,2), ss_alb_snw_drc(71,3),     &
            ss_alb_snw_drc(71,4), ss_alb_snw_drc(71,5)
         write (iulog,*) 'SNICAR: Mie single scatter albedos for diffuse ice, rds=100um: ',     &
            ss_alb_snw_dfs(71,1), ss_alb_snw_dfs(71,2), ss_alb_snw_dfs(71,3),     &
            ss_alb_snw_dfs(71,4), ss_alb_snw_dfs(71,5)
         IF (DO_SNO_OC) THEN
            write (iulog,*) 'SNICAR: Including OC aerosols from snow radiative transfer calculations'
         ELSE
            write (iulog,*) 'SNICAR: Excluding OC aerosols from snow radiative transfer calculations'
         ENDIF
      ENDIF
      !
#ifdef MODAL_AER
      IF (p_is_master) THEN
         ! unique dimensionality for modal aerosol optical properties
         write (iulog,*) 'SNICAR: Subset of Mie single scatter albedos for BC: ', &
            ss_alb_bc1(1,1), ss_alb_bc1(1,2), ss_alb_bc1(2,1), ss_alb_bc1(5,1), ss_alb_bc1(1,10), ss_alb_bc2(1,10)
         write (iulog,*) 'SNICAR: Subset of Mie mass extinction coefficients for BC: ', &
            ext_cff_mss_bc2(1,1), ext_cff_mss_bc2(1,2), ext_cff_mss_bc2(2,1), ext_cff_mss_bc2(5,1), ext_cff_mss_bc2(1,10),&
            ext_cff_mss_bc1(1,10)
         write (iulog,*) 'SNICAR: Subset of Mie asymmetry parameters for BC: ', &
            asm_prm_bc1(1,1), asm_prm_bc1(1,2), asm_prm_bc1(2,1), asm_prm_bc1(5,1), asm_prm_bc1(1,10), asm_prm_bc2(1,10)
         write (iulog,*) 'SNICAR: Subset of BC absorption enhancement factors: ', &
            bcenh(1,1,1), bcenh(1,2,1), bcenh(1,1,2), bcenh(2,1,1), bcenh(5,10,1), bcenh(5,1,8), bcenh(5,10,8)
      ENDIF
#else
      IF (p_is_master) THEN
         write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophillic BC: ', &
            ss_alb_bc1(1), ss_alb_bc1(2), ss_alb_bc1(3), ss_alb_bc1(4), ss_alb_bc1(5)
         write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophobic BC: ', &
            ss_alb_bc2(1), ss_alb_bc2(2), ss_alb_bc2(3), ss_alb_bc2(4), ss_alb_bc2(5)
      ENDIF
#endif

      IF (p_is_master) THEN
         IF (DO_SNO_OC) THEN
            write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophillic OC: ', &
               ss_alb_oc1(1), ss_alb_oc1(2), ss_alb_oc1(3), ss_alb_oc1(4), ss_alb_oc1(5)
            write (iulog,*) 'SNICAR: Mie single scatter albedos for hydrophobic OC: ', &
               ss_alb_oc2(1), ss_alb_oc2(2), ss_alb_oc2(3), ss_alb_oc2(4), ss_alb_oc2(5)
         ENDIF

         write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 1: ', &
            ss_alb_dst1(1), ss_alb_dst1(2), ss_alb_dst1(3), ss_alb_dst1(4), ss_alb_dst1(5)
         write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 2: ', &
            ss_alb_dst2(1), ss_alb_dst2(2), ss_alb_dst2(3), ss_alb_dst2(4), ss_alb_dst2(5)
         write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 3: ', &
            ss_alb_dst3(1), ss_alb_dst3(2), ss_alb_dst3(3), ss_alb_dst3(4), ss_alb_dst3(5)
         write (iulog,*) 'SNICAR: Mie single scatter albedos for dust species 4: ', &
            ss_alb_dst4(1), ss_alb_dst4(2), ss_alb_dst4(3), ss_alb_dst4(4), ss_alb_dst4(5)
         write(iulog,*)
      ENDIF

   END SUBROUTINE SnowOptics_init
   !-----------------------------------------------------------------------


   SUBROUTINE SnowAge_init( fsnowaging )

   USE MOD_NetCDFSerial

   IMPLICIT NONE

   character(len=256), intent(in) :: fsnowaging      ! snow aging parameters file name
   character(len= 32) :: subname = 'SnowAge_init'    ! SUBROUTINE name
      !
      ! Open snow aging (effective radius evolution) file:
      IF (p_is_master) THEN
         write(iulog,*) 'Attempting to read snow aging parameters .....'
         write(iulog,*) subname,trim(fsnowaging)
      ENDIF

      !
      ! SNOW aging parameters
      !
      CALL ncio_read_bcast_serial (fsnowaging, 'tau', snowage_tau)
      CALL ncio_read_bcast_serial (fsnowaging, 'kappa', snowage_kappa)
      CALL ncio_read_bcast_serial (fsnowaging, 'drdsdt0', snowage_drdt0)

      !
      IF (p_is_master) THEN
         write(iulog,*) 'Successfully read snow aging properties'
      ENDIF
      !
      ! print some diagnostics:
      IF (p_is_master) THEN
         write (iulog,*) 'SNICAR: snowage tau for T=263K, dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_tau(3,11,9)
         write (iulog,*) 'SNICAR: snowage kappa for T=263K, dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_kappa(3,11,9)
         write (iulog,*) 'SNICAR: snowage dr/dt_0 for T=263K, dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_drdt0(3,11,9)
      ENDIF

   END SUBROUTINE SnowAge_init
   !-----------------------------------------------------------------------


   real(r8) FUNCTION FreshSnowRadius (forc_t)
!-----------------------------------------------------------------------
! !DESCRIPTION:
!  Returns fresh snow grain radius, which is linearly dependent on temperature.
!  This is implemented to remedy an outstanding bias that SNICAR has in initial
!  grain size. See e.g. Sandells et al, 2017 for a discussion (10.5194/tc-11-229-2017).
!
!  Yang et al. (2017), 10.1016/j.jqsrt.2016.03.033
!  discusses grain size observations, which suggest a temperature dependence.
!
! !REVISION HISTORY:
!  Author: Leo VanKampenhout
!
!-----------------------------------------------------------------------
! !USES:
   USE MOD_Const_Physical, only: tfrz
   USE MOD_Aerosol, only: fresh_snw_rds_max

   ! !ARGUMENTS:
   real(r8), intent(in) :: forc_t         ! atmospheric temperature (Kelvin)
   !
   ! !LOCAL VARIABLES:
   !-----------------------------------------------------------------------
   real(r8), parameter :: tmin = tfrz - 30._r8       ! start of linear ramp
   real(r8), parameter :: tmax = tfrz - 0._r8        ! END of linear ramp
   real(r8), parameter :: gs_min = snw_rds_min       ! minimum value
   real(r8)            :: gs_max                     ! maximum value

      IF ( fresh_snw_rds_max <= snw_rds_min )THEN
         FreshSnowRadius = snw_rds_min
      ELSE
         gs_max = fresh_snw_rds_max

         IF (forc_t < tmin) THEN
            FreshSnowRadius = gs_min
         ELSEIF (forc_t > tmax) THEN
            FreshSnowRadius = gs_max
         ELSE
            FreshSnowRadius = (tmax-forc_t)/(tmax-tmin)*gs_min + &
                              (forc_t-tmin)/(tmax-tmin)*gs_max
         ENDIF
      ENDIF

   END FUNCTION FreshSnowRadius

END MODULE MOD_SnowSnicar
! ---------- EOP ------------
